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The  focus  of  this  work  is  the  derivation  and  development  of  a  class  of  models 
called  Generalized  Linear  FeedForvvard  (GLFF)  networks.  Included  is  a  description 
of  the  a  priori  assumptions  concerning  the  systems  under  study,  a  development  and 
formulation  of  the  models  (networks)  to  be  used,  and  a  definition  of  a  measurement 
criteria  used  to  judge  goodness-of-fit.  The  models  are  constructed  from  a  weighted 
sum  of  parametrically  dependent  kernels.  These  kernels  represent  basis  vectors  that 
span  a  vector  space  with  dimension  equal  to  the  order  of  the  network  realization. 
Optimization  of  the  model  parameters  is  discussed  and  a  new  approach  is  given;  one 
that  extends  the  conventional  squared  error  analysis  to  the  complex  domain. 

After  presentation  of  GLFF  model  theory  and  parametric  optimization 
methodologies,  network  realizations  are  constructed  to  address  the  system  modeling 
problem.      GLFF  networks  can  be  formulated  to  model  discrete-time,  continuous- 


time,  and  sampled-data  systems.  In  this  study,  it  is  shown  how  a  flight  motion 
simulator,  used  in  the  laboratory  test  and  evaluation  of  guided  weapon  systems,  can 
be  accurately  modeled  using  GLFF  networks  comprised  of  Gamma,  Laguerre,  and 
Legendre  kernels. 

Advantages  of  GLFF  networks  are  their  ability  to  accurately  represent  linear  time- 
invariant  systems  coupled  with  the  simplicity  of  the  procedures  employed  to  optimally 
determine  their  parameters.  The  network  structures  have  efficient  and  practical 
implementations  and  applications  are  offered  as  proof. 


VI 


CHAPTER  1 
PROBLEM  FORMULATION 


In  this  chapter  an  overview  of  the  entire  content  of  the  research  is  given, 
beginning  with  a  definition  of  the  problem  to  be  solved;  namely,  modeling  the 
impulse  and  step  responses  of  causal,  linear,  time-invariant  systems.  The 
system/model  configuration  is  illustrated  and  the  primary  measurement  objective 
stated.  The  model  architecture  is  briefly  described  and  prefaced  with  a  review  of  the 
major  contributors  leading  to  the  proposed  formulation.  Finally,  the  contribution  of 
this  research  is  summarized.  It  will  be  developed  throughout  the  remaining  chapters. 

It  will  be  shown  that  the  models  developed  arc  capable  of  approximating  systems 
with  real  and/or  complex  conjugate  pair  poles,  even  those  which  have  slowly 
decaying  impulse  responses  (large  time  constants)  and/or  oscillatory  modes  (light 
damping).  This  is  difficult  to  achieve  with  finite  impulse  response  (FIR)  structures. 
The  models  studied  are  considered  parametric  since  they  are  functions  of  parameter 
vectors  that  must  be  determined.  These  parameter  vectors  are  fairly  easy  to  optimize 
and  unlike  infinite  impulse  response  (IIR)  structures  stability  can  always  be 
guaranteed. 

The  networks  under  consideration  in  this  research  are  comprised  of  linear 
combinations  of  special  basis  fund  ions.  These  basis  functions  are  not  arbitrarily 
chosen  but  are  derived  from  assumptions  of  the  systems  being  modeled.     The 


networks  therefore  become  good  approximates  since  they  obtain  the  natural 
information  of  the  unknown  system.  This  information  resides  in  a  parameter  vector 
built  into  the  model  structure.  The  parameters  are  optimally  chosen  so  that  the 
models  approximate  the  system  in  a  least  squares  sense.  Hence,  the  proposed  theory 
is  an  implementation  of  parametric  least  squares  error  modeling  [Cad90]. 

The  modeling  methodology  can  also  be  described  using  geometric  concepts. 
Considering  signals  as  vectors  in  vector  spaces,  the  models  are  linear  combinations  of 
parametric  basis  vectors.  They  span  a  subspace  of  the  space  spanned  by  the  unknown 
system.  The  model  output  is  a  projection  of  the  system  vector  onto  the  model 
subspace.  This  space  has  an  orientation  that  is  parametrically  established.  The  error 
is  a  vector  that  lies  outside  the  model  subspace.  Projecting  orthogonally  the  system 
onto  the  model  space  minimizes  the  error. 

1.1  Problem  Definition 

The  objective  is  to  approximate  a  causal,  linear  time-invariant  (LTI)  system,  H. 
using  a  pre-defined  model,  M(p),  such  that  a  measure  of  the  quality  of  the 
approximation  is  minimized.  Models  developed  in  this  study  will  be  called 
generalized  linear  feedforward  (GLFF)  networks  and  will  be  defined  in  Chapter  3. 
The  basic  system  modeling  configuration  is  shown  in  Figure  1-1  below.  Arguments 
are  not  included  with  the  functions  {.v,  d,  v,  e,  H)  to  illustrate  that  the  formulation  will 
apply  in  both  a  continuous-time  setting  [x(t),  d(t),  y(t),  e(t),  H(s)}  and  a  sampled-data 
or  discrete-time  setting  [x(k),  d(k),  y(k),  e(k),  H(z)}.  Quite  often  in  this  thesis  we 
switch  back-and-forth  between  continuous-time  and  discrete-time  concepts.    This  is 


intentionally  done  to  prove  the  broad  applicability  of  this  work.  Most  of  the  theory 
however  will  be  derived  in  continuous-time.  It  would  have  been  just  as  appropriate  to 
develop  the  concepts  in  discrete-time. 
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Figure  1-1:  System  modeling  configuration. 


1.2  Measurement  Objective 

How  is  it  determined  whether  a  model  has  achieved  a  good  representation  o\'  the 
system  being  approximated?  For  this  a  measurement  o(  the  error  or  a  test  for 
goodness-of-fit  of  the  model  is  needed.  In  this  thesis  the  integrated  squared  error 
criteria  will  be  employed;  namely,  for  a  given  model  M(p),  the  parameter  set  p  is 
found  such  that  the  energy  in  the  error  signal  e  is  minimized.  This  is  one  of  the 
simplest  error  criteria  to  work  with  since  the  measurement  is  basically  quadratic  in 
form  and  its  derivatives  are  linear  in  certain  elements  of  the  parameter  set.  So  the 
objective  is  to  minimize 


or 


J=\  c'(t)cit=\  (d(t)-y(t))2dt     (continuous-time) 
J  =  ^V(£)  =  £(</(/;)-  v(/;)):   (discrete-time) 


(1.1 


*=0 


Jfc=0 


1.3  Model  Classifications 

To  provide  motivation  for  the  GLFF  networks  (a  complete  formulation  will  be 
given  in  Chapter  3),  consider  the  discrete-time  setting.  In  the  past  discrete-time 
system  models  have  been  largely  divided  into  two  classes:  FIR  and  IIR.  A  new  class 
of  models  (filters)  known  as  the  generalized  feedforward  (GFF)  filters  was 
popularized  with  the  development  of  the  Gamma  filter  [Pri93].  This  filter  has  a 
structure  similar  to  the  tap-delay  line  (FIR  filter)  with  the  delays  replaced  by  a  "more 
general"  transfer  function  known  as  the  Gamma  kernel.  Extending  this  development 
can  be  accomplished  by  considering  other  popular  kernels,  most  of  which  are 
orthogonal,  such  as  the  Laguerre,  Legendre,  Jacobi,  and  Kautz  [Lee67],  [Arm59], 
[Kau54].  The  extended  formulation  with  a  more  comprehensive  set  of  kernels  is 
called  the  GLFF  model  class  because  its  members  are  used  in  the  context  of  modeling 
linear  systems  and  they  obey  the  superposition  principle.  How  does  the  GLFF  model 
class  compare  to  FIR  and  IIR  models?  Figure  1-2  relates  these  three  classes  based 
upon  important  issues  that  are  usually  considered  when  deciding  which  modeling  tool 
to  pull  from  the  signal  processing  toolbox.  GLFF  models  can  be  said  to  fall 
somewhere  between  FIR  and  IIR  models  in  capability,  complexity,  stability  and 
computability. 

FIR  GLFF  IIR 

[complexity] 
low       <—       I  >       — >        high 

[  capability  J 

f      stabilitv      1 
trivial     <—     <  >    — >     difficult 

[computability  I 

Figure  1-2:  FIR/GLFF/IIR  comparison. 


This  thesis  exploits  the  strengths  of  GLFF  models.  They  offer  an  alternative  to 
the  two  extremes  in  the  modeling  process.  Their  advantages  are  that  they  have  a 
reduced  complexity  compared  to  IIR  models  while  offering  greater  capability  over 
FIR  models.  Unlike  IIR  modeling  techniques,  stability  can  be  easily  guaranteed  and 
their  optimal  or  near-optimal  solutions  can  be  computed  with  much  less  difficulty. 
The  most  difficult  task  is  in  selection  of  the  time  scale  vector  A.  This  is  one  element 
of  the  parameter  set  p.  The  complete  parameter  set  will  be  defined  in  Chapter  3. 
Optimization  of  its  elements  comprises  all  of  Chapter  4. 

1.4  Historical  Review 

Much  of  the  theory  on  which  the  GLFF  model  class  is  founded  began  with  the 
work  of  Lee  [Lee33],  [Lee67].  Lee  showed  how  to  use  orthogonal  functions  to 
synthesize  transient  networks  and  how  this  concept  is  employed  to  build  electrical 
devices  that  can  realize  physical  systems.  Another  major  contributor  is  Kautz 
[Kau54],  He  provided  a  convenient  way  to  develop  the  time  domain  response  of  a 
system  using  a  most  general  orthogonal  signal  set.  Young  and  Huggins  extended 
Kautz'  procedure  to  discrete-time  [You62a]  as  did  Broome  [Bro65].  Nurges  and 
Yaaksoo  [Nur81]  formulated  state  space  equations  that  represent  discrete-time  linear 
systems  using  the  Laguerre  model.  Roberts  and  Mullis  provide  a  generalized  state 
variable  formulation  of  orthogonal  all-pass  network  structures  [Rob87].  These 
structures  keep  data  paths  local.  The  fundamental  processing  elements  are  alike 
except  perhaps  for  local  scaling  parameters.  This  simplifies  the  design  and  allows  for 
distributing  the  computations  in  a  connected  array.    These  orthogonal  structures  also 


have  the  advantage  that  higher  order  models  can  be  constructed  by  simply  adding  new 
blocks  or  terms  without  effecting  or  forcing  a  recalculation  of  the  existing  lower  order 
terms.  The  concept  of  a  generalized  feedforward  filter  was  proposed  by  Principe 
[Pri93]  with  the  Gamma  model.  This  idea  of  creating  a  more  general  transversal  filter 
using  the  earlier  work  with  orthogonal  signal  sets  leads  to  the  modeling  methodology 
here  described;  namely,  the  generalized  linear  feedforward  networks. 

What  is  the  basic  objective?  Given  a  desired  impulse  response,  h(t),  of  a  causal 
LTI  system  along  with  a  prescribed  tolerance  on  the  allowable  error,  find  an 
approximating  impulse  response,  h(t),  that  is  simple  in  structure  and  easy  to  compute. 
Let  the  system  and  model  transfer  functions  be  H(s)  and  H(s),  respectively. 
Construct  H(s)  with  the  form 

m 

H(s)=m=AWs~z,)  =±-A_  (1.2) 

1=1 
Assuming  no  coincident  poles  the  impulse  response  is 

1=1 

where  Re(/?,)  <  0,  V  p=1,  ...,  n.  It  is  desired  that  h(t)  and  h(t)  be  as  close  as  possible 
under  prescribed  tolerances. 

There  are  two  well-known  approaches  to  performing  this  approximation.  The  first 
is  to  transform  h(t)  into  H(s)  and  then  approximate  H(s)  in  the  .v-domain  by  using  a 
ratio  of  polynomials.  This  approach  is  generally  less  difficult  and  several  techniques 
for  accomplishing  this  approximation  have  been  presented  in  the  literature.    These 


include  Prony's  method,  Pade  approximants.  continued-fraction  expansions,  and  even 
analog  filter  design  [Su71],  [Bak81],  [Opp89].  The  problem  with  this  approach  is 
that  the  error  in  the  time  domain  is  difficult  to  control.  The  second  approach  is  to 
approximate  h{t)  with  a  finite  sum  of  exponentially  damped  sinusoids  and  damped 
exponentials  (as  described  above)  and  then  transform  the  approximation  to  find  H(s). 
Here  the  modeling  error  may  be  controlled  directly  in  the  time  domain.  One  way  of 
accomplishing  this  is  to  implement  the  procedure  described  by  Kautz  [Kau54]  which 
involves  approximating  h(t)  using  a  weighted  sum  of  special  basis  functions.  These 
basis  functions  are  themselves  weighted  sums  of  exponentials  carefully  constructed  so 
that  they  are  orthonormal.  This  provides  for  great  simplification  when  searching  for 
the  optimal  parameter  set  p  and  gives  direct  control  over  the  associated  error. 

1.5  Modeling  Methodology 

In  this  thesis  the  method  used  to  implement  GLFF  network  approximations  of 
linear  systems  is  based  on  the  work  done  by  Kautz.  This  is  earned  out  by  considering 
a  construction  of  /;(/)  using  an  infinite  sum  of  weighted  time  domain  basis  functions 
{<pn(t,X  )}  .  Since  the  infinite  sum  must  be  truncated  to  a  finite  number  of  terms  an 
approximation  results 

//(/)  =  //(/)=  £u;,0„(/,/l)  (1.4) 

n=0 

The  weights  vv„,  w\ wN  and  basis  functions  {0„(t,X)\      are  chosen  with  the 

following  goals  in  mind: 

1.    For  a  fixed  number  of  terms,  N,  make  the  convergence  as  rapid  as  possible. 
This  is  controlled  using  a  time  scale  vector  X. 


2.  Make  the  optimal  weights  independent  of  the  number  of  taps,  N+\.  This  is 
accomplished  by  using  orthogonal  basis  functions. 

3.  For  a  fixed  number  of  taps,  N+l,  minimize  the  energy  in  the  error  signal.  This 
requires  determining  optimal  values  for  all  of  the  elements  of  the  parameter 
setp  in  such  a  way  that  J  (the  integrated  squared  error)  is  minimized. 

The  above  procedure  is  applicable  when  modeling  the  impulse  response  of  a 

desired  system.    It  may  also  be  of  interest  to  approximate  the  step  response.    In  fact 

the  applications  given   in  this  thesis  involve  finding  system  models  using  step 

response  data.    If  the  network  input  is  a  step  function,  this  method  is  employed  by 

approximating  not  the  desired  system  output  but  its  derivative  or  by  approximating 

the  desired  output  and  taking  the  derivative  of  the  result.    This  will  give  an  analytic 

expression  for  the  impulse  response.  When  only  a  numerical  or  tabulated  form  of  the 

impulse  response  is  desired,  use  the  step  response  data  to  generate  a  system  model. 

Once  the  model  is  established,  use  an  impulse  (or  white  noise)  input  to  generate  the 

impulse  response  data. 

1.6  Research  Contribution 

This  study  extends  the  above  concepts  and  contributes  to  the  system  modeling 
theory  in  several  ways. 

First,  a  unified  system  modeling  framework  centered  around  the  GLFF  model 
class  is  provided.  As  will  be  seen  in  Chapter  3  this  formulation  not  only  includes 
models  that  employ  orthogonal  signal  sets  like  the  Laguerre  and  Kautz  but  non- 
orthogonal  ones  as  well,  such  as  the  Gamma.  In  addition  to  these  generalized  models, 
the  basic  framework  can  also  be  used  to  represent  FIR  and  IIR  structures. 


Optimal  parameter  selection  is  a  big  part  of  the  modeling  procedure.  The  most 
difficult  task  is  locating  or  approximating  the  optimal  time  scale  vector.  A  new 
approach  that  extends  traditional  squared  error  analysis  to  the  complex  domain  is 
developed.  It  will  be  shown  how  the  information  in  this  new  space  can  be  used  to 
locate  the  optimal  time  scale  parameters. 

Finally  an  application  of  the  theory  is  given  by  modeling  the  step  response  of  a 
Carco  5-axis  flight  motion  simulator  (FMS).  This  system  is  used  in  hardware-in-the- 
loop  (HITL)  testing  of  guided  weapon  systems.  GLFF  modeling  of  other  systems  and 
subsystems  used  in  HITL  testing  is  also  considered.  It  is  demonstrated  how  this  class 
of  models  can  be  employed  to  accurately  represent  the  dynamics  of  physical  systems 
and  add  value  to  the  task  of  weapon  system  test  and  evaluation. 


CHAPTER  2 
SYSTEM  ASSUMPTIONS 


In  the  previous  chapter,  a  brief  discussion  was  given  concerning  the  assumed 
systems  to  be  modeled.  A  more  complete  analysis  of  all  the  system  assumptions  is 
now  offered.  Mathematical  and  physical  arguments  are  used  to  validate  their 
necessity.  The  basic  assumptions  are  that  input/output  signals  have  finite  energy  and 
systems  are  stable,  causal,  linear,  time-invariant,  and  realizable.  After  stating  the 
assumptions,  the  technique  employed  to  synthesis  these  systems  by  simulation  is 
given.  Higher  order  systems  will  be  simulated  by  cascading  together  1st  and  2nd  order 
subsystems.  This  is  done  in  a  continuous-time  framework.  The  transformation  used 
when  a  digital  system  is  desired  is  also  discussed.  Finally  some  simulation  examples 
of  the  kinds  of  real  systems  to  be  modeled  are  offered. 

2.1  Linear,  Time-Invariant,  Causal  Systems 

There  are  many  ways  to  classify  systems;  linear  vs.  nonlinear,  time-invariant  vs. 
time  varying,  causal  vs.  non-causal,  stochastic  vs.  deterministic,  continuous  vs. 
discrete,  etc.  One  important  step  in  the  study  and  analysis  of  systems  theory  is  to 
derive  a  structural  or  mathematical  description  of  processes  and  systems.  This 
assumed  structure  is  essential  if  plausible  analytical  models  are  sought  to  approximate 
them.  One  could  argue  that  if  we  had  a  mathematical  description  of  the  process  why 
would   we   need   additional    models.      The   primary  reason   for  the   model    is   to 
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approximate  input/output  behavior  of  the  process.  The  assumed  mathematical 
description  of  the  process  could  be  used  to  give  very  specific  information  about  its 
internal  composition  (time  constants,  damping  ratios,  etc.);  however,  a  component 
level  specification  may  be  unavailable.  In  fact  even  if  these  system  parameters  are 
known  they  may  not  be  operating  within  their  stated  specification.  In  this  research  an 
understanding  of  the  internal  composition  is  not  the  objective,  only  a  model  of  the 
input/output  behavior  when  viewing  the  system  as  a  black  box.  The  theory  presented 
in  this  section  does  however  give  the  foundation  used  to  formulate  and  validate  this 
new  class  of  models  and  serves  to  illustrate  the  derivation  of  these  concepts. 

Many  systems  studied  in  engineering  can  be  well  represented  by  mathematical 
relations  that  describe  their  physical  makeup  or  behavior.  Quite  often  systems  are 
very  complex  and  a  realistic  characterization  may  require  nonlinear  and/or  time- 
varying  equations  that  are  difficult  to  solve.  For  simplicity,  practicality,  and  in  order 
to  establish  an  applicable  set  of  tools  which  can  be  used  in  the  design  and  analysis  of 
systems,  approximations  and  assumptions  are  made  about  their  physical  nature.  This 
is  done  primarily  to  allow  for  the  use  of  linear  systems  theory.  This  reasoning  is  often 
justified  in  one  of  two  ways: 

1 .  The  system  is  assumed  to  be  linear  by  nature  or  it  at  least  operates  in  a  linear 
region  of  its  capabilities  so  that  the  linearity  conditions  are  satisfied. 

2.  The  system  is  assumed  to  be  nonlinear  by  nature  and  in  order  to  apply  linear 
systems  theory  we  conceptualize  a  linearization  around  a  nominal  operating 
point.  With  this  approach  precautions  must  be  taken  to  ensure  the  domain  of 
the  linear  assumption  is  not  exceeded. 

Why  is  the  use  of  linear  systems  theory  desirable?    In  addition  to  the  availability 

of  a  well-studied  set  of  tools,  most  often  a  description  of  a  system  via  the  impulse 
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response  and  transfer  function  is  desired.  These  functional  system  descriptions  are 
established  by  employing  the  theory  of  linear  systems  [DeC89]. 

Consider  a  linear  time-invariant  (LTI)  system  having  input  x{t)  and  output  d(t). 
The  impulse  response,  h(t),  is  defined  as  the  output  when  the  input  is  a  unit  impulse 
function  d\t).  The  transfer  function,  H(s),  is  defined  as  the  Laplace  transform  of  h(t) 
with  zero  initial  conditions.  This  is  equivalent  to  assuming  h(t)  =  0  for  t  <  0  (causal). 

The  above  definition  of  the  transfer  function  is  in  terms  of  the  impulse  response. 
However  a  description  of  the  transfer  function  taking  a  different,  more  mathematical 
approach  is  also  available.  In  practice,  a  large  class  of  signals  and  systems  in 
engineering  can  be  formulated  through  the  use  of  differential  equations.  The  LTI, 
causal  assumptions  equivalently  translate  into  the  requirement  that  systems  be 
representable  by  linear,  constant  coefficient,  ordinary  differential  equations.  Namely, 
an  n    order  system  can  be  written  as 

dw(t)  +  an_ld{n-1\t)  +  -  +  a1dll\t)  +  aQd(t)  = 

bmx{m)(t)  +  bm_]x{m-l)(t)  +  ---  +  bi.x{U(t)  +  b0.x(t) 

where  a^  a\,  ...,  an.\,  bo,  b\,  ...,  bm  are  all  real  constants  with  n  >  m.  Systems 
described  by  ordinary  differential  equations  of  this  form  are  said  to  contain  lumped 
parameters  rather  than  distributed  parameters  [Kuo87].  Partial  differential  equations 
such  as  the  elliptic  equations,  hyperbolic  equations,  and  parabolic  equations  describe 
phenomena  such  as  fluid  flow,  wave  propagation,  and  heat  conduction  all  of  which 
are  considered  to  be  distributed  parameter  systems  [Pow79]. 
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To  obtain  the  transfer  function  of  an  LTI  system  described  by  the  differential 

equation   above  assume  zero  initial  conditions   (causality)  and  take  the  Laplace 

transform 

II{s)=Ws)  =  l>ms'»+bm_lsm-x+-+bls  +  b6  (22) 

X(s)       sn  +  an_ts"'1  +■  •  -+a,  s  +  a0 

Rewriting  this  in  factored  form  gives 

m 

IK*-*) 

ff(s)  =  A-f (2.3) 

1  =  1 

and  in  partial  fraction  form  (assuming  no  repeated  poles) 

!  i  {s-Pi) 

To  obtain  the  impulse  response,  take  the  inverse  Laplace  transform.    Assuming  no 
coincident  poles  this  yields 


/»(')  =  £vw  (2.5) 


1=1 


So  the  LTI,  causal  assumption  leads  to  the  following  result:  LTI  causal  systems 
can  be  represented  by  constant  coefficient  ordinary  differential  equations  whose 
solutions  (impulse  responses)  are  comprised  of  a  sum  of  exponential  functions  that 
are  real  (when  p,  is  real)  and/or  complex  (when  p,  is  complex).  If  some  of  the  poles 
are  repeated  then  there  are  terms  with  multiplicative  factors  t,  r,  r\  ...,  r  for  a  k 
order  pole.  This  representation  is  the  motivation  for  the  modeling  structure 
considered  in  this  thesis. 
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2.2  Realizubility 

It  was  shown  that  the  LTI,  causal  system  assumption  leads  to  a  transfer  function 
representation  of  the  form 

m 

n('-«) 

S(J)=A-? (2.6) 

An  answer  is  now  given  to  the  following  question  "What  are  the  constraints  on 
this  form  that  ensures  it  is  a  description  of  a  physically  realizable  system?"  The 
realizability  requirement  of  H{s)  is  satisfied  with  the  following  constraints: 

1 .  The  constant  A  must  be  a  real. 

2.  The  number  of  poles,  /?,  must  be  greater  than  the  number  of  zeros,  m. 

3.  Complex  poles  and  zeros  must  occur  in  conjugate  form.  That  is  if  p,  is 
complex  then  there  must  be  a/?,,  for  some  j,  such  that  p]  =  p*  where  *  denotes 
complex  conjugation. 

When  H(s)  is  realizable,  it  is  easy  to  show  that  the  impulse  response  is  represented 

by  a  weighted  sum  of  exponentials  and  exponentially  decaying  sinusoids. 

2.3  L~  Signals  and  Systems 

The  final  constraint  imposed  is  that  the  signals  must  belong  to  the  vector  space 
L"(IR),  the  square-integrable  functions.  Here  integrable  is  meant  in  the  Lebesgue 
sense  [Deb90].  Assume  the  signal  xe  L2([R)  and.v:  R  \-^  R  then 


]x2(t)dt«*>  (2.7) 


Hence  it  is  said  that  x  has  finite  energy.  Let  h(t)  be  the  impulse  response  of  a  system. 
Since  h(t)  =  0  for  t  <  0,  and  he  L2(R)  then 
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\h\t)dt  <oo  (2.8) 


In  this  case  it  is  said  that  h  is  a  stable  system.  So  by  operating  in  L"(K)  a  restriction  is 
imposed  to  work  only  with  finite  energy  signals  and  stable  systems. 

2.4  System  Synthesis 

Recall  the  system  realizability  requirement  results  in  an  impulse  response 
representation  that  is  a  weighted  sum  of  exponentials  and  exponentially  decaying 
sinusoids.  From  the  transfer  function  perspective  then  a  system  can  be  thought  of  as 
consisting  of  a  cascade  of  Is1  and  2nd  order  subsystems.  The  most  general  forms  of 
these  subsystem  structures  are 


H,W  =  ^M- 


X(s)      s  +  aw 
HJs)  = 


D2(s)  _      b2]s  +  b2() 


(2.9) 


X(s)      r  +  a21s  +  a2i) 
resulting  in  differential  equations 

dl(t)+ai0dl(t)  =  bl0x(t) 

d2(t)  +  a2]d2(t)  +  a2„d{t)  =  b2lx(t)  +  b2i)x(t) 


(2.10) 


During  this  thesis,  examples  will  be  given  whereby  simulated  systems  are  created. 
For  such  examples,  the  synthesized  systems  will  consist  of  these  P1  and  2nd  order 
subsystem  blocks.  Most  often  slightly  less  general  forms  will  be  used  to  keep  the 
synthesis  more  in  line  with  engineering  (rather  than  mathematical)  concepts.  These 
special  forms  convey  the  engineering  information  of  usual  interest:  time  constants, 
(lamping  ratios,  and  undamped  natural  frequencies.  These  subsystems  have  the 
following  transfer  functions 
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ffiW...S«__L.__2L 

x(s)    ts+1    s+yT 


(2.11) 


X(s)      s2  +  2gct)ns  +  (0] \ 
and  corresponding  differential  equations 

Td1(t)  +  dl(t)  =  x(t) 
d2(t)  +  2^a)nd2(t)  +  Q)2nd(t)  =  0)2nx(t) 

where 


(2.12) 


r  =  time  constant 

£   =  damping  ratio 

(On  -  undamped  natural  frequency 

2.5  Continuous/Discrete-Time  Transformation 

The  goal  of  this  thesis  is  to  define  a  model  class  useful  for  approximating  discrete- 
time,  continuous-time,  and  sampled-data  systems.  The  above  system  formulation  was 
carried  out  in  the  continuous-time  domain.  The  same  principles  however  could  have 
been  derived  in  the  discrete-time  domain.  All  the  concepts  discussed  have  equivalent 
discrete-time  representations.  Sometimes  a  continuous-time  system  in  analytic  form 
is  given  when  a  discrete-time  representation  is  desired.  This  can  be  resolved  by 
sampling.  In  this  case,  the  discrete-time  system  is  called  a  sampled-data  system. 
Hence  a  sampling  time,  T,  will  necessarily  be  specified. 

Having  a  continuous-time  transfer  function,  there  are  several  methods  used  to  map 
to  the  discrete-time  domain.  The  most  popular  are  the  bilinear,  impulse-  or  step- 
invariant,  and  matched-Z  transform  [Str88].  There  are  tradeoffs  among  these 
mappings  and  some  work  well  where  others  do  not.  The  most  widely  used  is  the 
bilinear  transform.  The  transformation  (or  mapping)  is  carried  out  as  follows 
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H(z)=H(s)[__2S-1  (2.13) 


T  :  +  \ 


where  T  is  the  sampling  interval.  This  is  also  known  as  Tustin's  method  in  control 
systems  circles  and  is  derived  by  considering  the  problem  of  numerical  integration  of 
differential  equations  [Fra90].  A  derivation  of  the  bilinear  transform  is  as  follows. 
Beginning  with  the  transfer  function  H(s),  determine  its  differential  equation 
representation.  A  difference  equation  is  then  derived  whose  solution  matches  the 
solution  of  the  differential  equation  at  the  sample  points.  The  specific  transformation 
depends  on  the  method  with  which  the  differential  equation  solution  is  computed. 
Consider  the  1st  order  system 

W(S)=A(£)=_JL_  (214) 

X(s)     s  +  yT 

or  associated  differential  equation 

4(t)+±4(')=ix(0  (2.15) 

The  solution,  d\(t),  can  be  written  in  integral  form  as 

dl(t)  =  ±j[-dl(v)  +  x(v)]dv  (2.16) 

0 

For  the  sampling  interval  of  time  T  this  can  be  rewritten  as 

nT-T  nT 

4(n)  =  4(fi7)  =  ±  j[-d,(v)  +  x(v)]dv+±  \[-dx{v)  +  x(v)]dv 

0  nT-T 

nT 

=  4(n-l)+-j:  \[-dl(v)  +  x(v)]tlu 


(2.17) 


nT-T 


The  method  in  which  the  integral  is  computed  over  the  region  [nT-T,  nT] 
determines  the  discrete-time  representation,  cl(n),  of  d{t)  and  hence  the  transformation 
used  to  map  H(s)  to  H(z).  Two  of  the  most  common  integration  schemes  are  Eider's 
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rule  (backward  difference  method)  and  the  trapezoidal  rule.   For  these  two  cases  the 


following  mappings  result. 


Table  2-1:  Continuous-time  to  discrete-time  mappings. 


Integration  Rule 

Mapping  H(s)  h>  H(z) 

Euler's 

A    t    J 

Trapezoidal 

rU+iJ 

So  the  bilinear  transform  is  derived  by  computing  the  integral  over  [nT-T,  nT\ 
using  the  trapezoidal  rule. 

In  Steiglitz's  work  [Ste65]  a  more  theoretical  treatment  and  significance  of  the 
bilinear  transform  is  developed.  In  fact,  it  is  shown  that  the  bilinear  transform  is  an 
isomorphism  connecting  the  continuous-time  vector  space  L  (IRL)  and  the  discrete-time 

(or  digital)  vector  space  of  double-ended  sequences  of  complex  numbers  {^„}°°=  x  that 

are  square  summable.  This  connection  proves  that  the  signal  processing  theories  with 
respect  to  LTI  causal  signals  and  systems  are  identical  in  both  the  continuous-time 
domain  and  discrete-time  domain.  In  this  thesis,  when  a  continuous-to-discrete 
transformation  is  performed  the  bilinear  transform  will  always  be  used. 

2.6  Examples 

Several  examples  of  the  1st  and  2nd  order  systems  discussed  in  Chapter  2.4  are 
shown  on  the  next  two  pages  in  Figures  2-1  to  2-4.  These  figures  illustrate  the 
impulse  and  step  response  behavior  that  must  be  approximated  by  the  class  of  models 
proposed  in  this  thesis. 


19 


1  8 

1 !                  I                  1                   1 

1                1                1                1 

16 

- 

14 

\    t=0  5 

- 

12 

-      \ 

- 

1 

\r=l  \ 

- 

0.8 

- 

0.6 

T=2           \. 

- 

04 

L_J^^~""~~---V>. 

- 

02 

1                 1                 1              1      ■ 1 

- 

n 

_     "1             — 1 1 — I           — 

05  1  15  2  25  3  3  5  4  45  5 

I  (sec) 

Figure  2-1:  Is1  order  impulse  responses,  7^=0.5,  1,  2,  3. 
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Figure  2-2:  1st  order  step  responses,  £=0.5,  1,  2,  3. 
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>nd 


Figure  2-3:  2    order  impulse  responses  with  ftj,=l. 
a)  £=0.4,  under  damped;  b)  £=1.0,  critically  damped;  c)  £=1.5,  over  damped. 
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Figure  2-4:  2     order  step  responses  with  <^,=  1. 
a)  £=0.4,  under  damped;  b)  £=1.0,  critically  damped;  c)  £=1.5,  over  damped. 
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2.7  Conclusion  (System  Assumptions) 

In  this  chapter  the  assumptions  made  about  the  systems  to  be  modeled  have  been 
stated;  causality,  linearity,  time-invariance,  realizability.  This  discussion  was  to 
provide  the  mathematical  foundation  and  engineering  principles  that  motivate  the 
development  of  the  system  models  to  be  introduced.  For  the  purposes  of  simulating 
real  systems,  1st  and  2nd  order  subsystems  are  used  as  the  basic  building  blocks.  These 
subsystems  can  be  cascaded  together  to  create  high  order  systems.  Methods  were 
discussed  that  allow  for  transforming  continuous-time  system  descriptions,  H(s),  into 
discrete-time  descriptions,  H{z). 

A  good  model  is  one  with  a  structure  that  can  capture  the  dynamics  of  a  system 
(time  constants,  damping  coefficients,  natural  frequencies)  without  explicit 
knowledge  of  its  parameters.  For  systems  with  large  time  constants  (slowly  decaying 
impulse  response)  and/or  small  damping  coefficients  (oscillatory  modes)  a  good 
approximation  is  more  difficult  to  achieve.  The  models  presented  in  this  thesis  are 
good  approximates  of  systems  even  with  these  characteristics. 


CHAPTER  3 
MODEL  FORMULATION 


In  this  chapter  our  first  major  contribution  to  the  problem  of  system  modeling  is 
presented;  namely,  a  comprehensive  framework  given  by  the  GLFF  model  class. 
Although  many  of  the  components  have  been  around  for  quite  some  time  it  is  the 
unifying  approach  taken  here  that  offers  new  insight.  This  new  formulation 
encompasses  not  only  GLFF  networks  but  FIR  and  IIR  models  (for  the  discrete-time 
case)  as  well. 

3.1  General  Model  Form  M(p) 

Recall  the  general  system  modeling  configuration,  shown  in  Figure  3-1  below 
where  H  is  the  system  to  be  approximated  using  the  model  M(p).  The  specific  form 
of  the  model  M{p)  is  now  discussed. 


d 

H 

v 

d) 

V 

€ 

M(p) 

y 

Figure  3-1:  System  modeling  configuration. 
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p  is  defined  as  a  parameter  set.  Its  elements  are  vectors  called  the  parameter 
vectors  of  the  model  M(p).  The  parameter  set  contains  three  elements  (vectors), 

p  =  {N,w,A] 

where  (3.1) 

N     =    dimension  vector 

vv     =     weight  vector 

A     =     time  scale  vector 

The  model  can  be  written  in  the  form 

M(p)=  M({N,W,A})  =  f(N,w,Hn(.tA))  (3.2) 

where  //„(•,  A)  is  called  the  model  kernel  at  stage  n,  0  <  n  <  N„  and  N,  is  a  component 
of  the  vector  N.  So  M(p)  is  a  function,  /,  of  the  parameter  vectors  and  the  model 
kernel.  The  kernel  argument  •  is  used  to  imply  that  the  domain  of  operation  (s  or  z) 
is  universal.  When  solving  specific  problems  the  notation  H„(s,A)  is  used  in  the 
continuous-time  case  and  H„(z,A)  in  the  discrete-time  case. 

3.2  Discrete-Time  Models 

Consider  the  definitions  above  in  the  discrete-time  domain.  This  unifying 
representation  can  be  used  to  illustrate  many  specific  model  structures.  In  this  section 
it  is  shown  how  to  describe  the  FIR,  IIR,  and  GLFF  model  classes  using  this 
formulation. 

3.2.1  FIR  Model 

In  the  case  of  the  FIR  models 

M(p)=M({N,n<,A})=fjn<iH„(z,A)  (3.3) 

where 
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N  =  (Nl) 

w  =  (w0,w1,...,wNf 


A  =  0  (3-4) 

Hn(zM  =  Hn(z)  =  flHi(z) 


i=0 


The  functions  H^z)  are  given  by 


/  =  0 
/>0 


HM  =  \    _,      .     „  (3.5) 


3.2.2  IIR  Model 


In  the  case  of  the  IIR  models 


2>A(^) 

M(p)  =  M({7V,  w,  A  })  =  -**r (3.6) 


n=l 


where 


N  =  (Nl,N2f 

=  (ava2,...,aNi,b0,b],...,bNt) 


w 

A  =  0 


Hn(z,A)=Hn(z)  =  Y\H,(z) 


i=0 


Just  as  in  the  FIR  case  the  functions  H((z )  are  given  by 


(3.7) 


_  f  1      /=0 


3.2.3  GLFF  Model 


In  the  case  of  the  GLFF  models 


M(p)  =  M({N,w,A  })  =  ]£  H'„//„(z,i )  (3.9) 

where 
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N  =  (N]) 

w  =  (w0,wr...,wN)T 

X  =  (AVAZ,...,AL)T  <3-'°> 

The  functions  HlX(z,A )  and  Hn2(z,A)  could  be  any  transfer  function.   In  this  thesis  a 
restriction  is  made  to  consider  special  cases  that 

1.  are  founded  on  the  system  model  assumptions  described  in  Chapter  2. 

2.  have  kernels  that  are  (typically)  orthogonal. 

3.  have  kernels  that  are  cascaded  lowpass  and  allpass  subsystems. 

4.  will  allow  for  modeling  the  unknown  system  with  real  and/or  complex  poles 
via  the  time  scale  vector  A. 

5.  have  localized  feedback. 

6.  guarantee  stability. 

3.2.4  FIR,  IIR,  GLFF  Model  Block  Diagrams 

To  compare  the  structures  of  these  models  examine  their  block  diagrams  in 
Figures  3-2,  3-3,  and  3-4.  As  can  be  seen  from  the  FIR  model,  its  lack  of  any 
feedback  component  limits  its  applicability.  The  IIR  model  offers  the  greatest 
capability  but  at  a  significant  cost  and  loss  of  stability  control.  The  GLFF  model 
offers  a  practical  alternative  that  rivals  the  IIR  model  in  performance  while 
maintaining  simple  computational  procedures  similar  to  the  FIR  case. 
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Figure  3-2:  FIR  model  block  diagram. 
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Figure  3-3:  IIR  model  block  diagram. 
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Figure  3-4:  GLFF  model  block  diagram. 
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3.3  GLFF  Model  Kernels 

The  GLFF  kernels  can  be  comprised  of  almost  any  transfer  function.  Some  of  the 
most  popular  discrete-time  and  continuous-time  kernels  are  given  in  Tables  3-1  and  3- 
3,  respectively.  These  tables  include  each  kernel's  common  name,  symbol,  and  nlh  tap 
transfer  function,  //„(•,  X).  Most  GLFF  networks  are  derived  from  orthogonal  signal 
sets  (see  Appendix  A).  This  orthogonality  of  a  kernel  is  also  indicated  in  the  tables. 
An  important  part  of  a  GLFF  network  is  the  time  scale  parameter  vector  A.  It  allows 
for  local  feedback  and  is  responsible  for  the  popularity  and  usefulness  of  this  entire 
class  of  models.  The  dimension  of  the  time  scale  vector,  dim(A),  determines  the 
number  of  degrees-of-freedom  a  kernel  offers.  dim(A)  is  also  given  in  the  tables. 
When  fitting  a  model  to  an  unknown  system,  the  parameter  set  must  be  optimally 
selected  to  achieve  the  minimum  integrated  squared  error. 

3.3.1  Popular  Discrete-Time  GLFF  Kernels 

The  most  popular  discrete-time  GLFF  kernels  are  the  delay,  Gamma,  Laguerre, 
Legendre,  and  Kautz.  They  are  listed  in  Table  3-1,  and  restrictions  on  the  time  scale 
parameters  are  given  in  Table  3-2.  Some  of  the  earliest  discussions  of  these  kernels 
arc  found  in  the  work  by  Young  [You62a],  Perez  [Per88],  and  Principe  [Pri93]. 

3.3.2  Popular  Continuous-Time  GLFF  Kernels 

The  most  popular  continuous-time  GLFF  kernels  are  the  Gamma,  Laguerre, 
Legendre,  Jacobi,  and  Kautz.  They  are  listed  in  Table  3-3,  and  restrictions  on  the  time 
scale  parameters  are  given  in  Table  3-4.  Some  of  first  to  describe  these  kernels  are 
Lee  [Lee33],  Kautz  [Kau54],  Armstrong  Ar[Arm59],  and  Principe  fPri93]. 
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Table  3-1:  Discrete-time  GLFF  kernels. 


Kernel 

nlh  tap 
symbol 

Hj2(z,A) 

Hkl{z,X) 

ff„(z,A) 

ortho- 
gonal 

dim(A) 

Delay 

H„(z) 

1 

z 

-/I 

z 

no 

0 

Gamma 

Gn{zj) 

1 

Az~] 

r    iz-1    ] 

« 

no 

1 

l-(i-A)z'1 

Vl-(l-^)z"'j 

Laguerre 

UzJ) 

Vi-/t2 

i-Az'1 

zx-A 

\-Az~1 

Vi-/i2 

1-az"1 

fl 

yes 

1 

Legendre 

P„(z,A) 

1 

z~]-Ak 
1-yLV1 

1          j-r.          A 

yes 

>  1 

\-A'zx 

Kautz 

Kn(z,A) 

^ 

*■> 

2 

V'-W2 

yes 

>  1 

1-/1  jf' 

i- 

-4.*"1 

NO    1-A^ 

Table  3-2:  Restrictions  on  discrete-time  time  scale  parameters. 


Kernel 

Valid  time  scale 

Delay 

none 

Gamma 

A  real,   0  <  A  <  2 

Laguerre 

A  real,    A  <  1 

Legendre 

Ak  real,  Ak  =Ak,0<A<\ 

Kautz 

Ak  complex,    Ak  <  1 
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Table  3-3:  Continuous-time  GLFF  kernels. 


Kernel 

na  tap 

symbol 

Hj2(s,A) 

Hkl(s,A) 

Hn(s,A) 

ortho- 
gonal 

dimU) 

Gamma 

Gn{s.A) 

1 

A 
s  +  A 

(    A    > 

Ks  +  A  , 

n 

no 

1 

Laguerre 

Ln(s,A) 

V2l 
s  +  A 

s-A 

s  +  A 

V2I 
s  +  A 

fs-A\' 
^s  +  AJ 

yes 

1 

Legendre 

P„(s,A) 

s  +  Aj 

s  -  Ak 

s  +  A, 

s  +  AnUs  +  A, 

yes 

>  1 

Jacobi 

Jn(sM 

s  +  Aj 

s-Ak 
s  +  Ak 

^t\s-A, 
s  +  Ani}s  +  Ai 

yes 

>  1 

Kautz 

Kn(s,A) 

ffiH&fiJp 

(s-Ak)(s-Ak) 

J2U4^ 

Pp-A,ls-t) 

yes 

>  1 

(s+AjXs+Aj) 

(s+Ak)(s+Ak) 

(s+AnXs+An) 

!£ 

MK.v+^) 

Table  3-4:  Restrictions  on  continuous-time  time  scale  parameters. 


Kernel 

Valid  time  scale 

Gamma 

A  real,  A>0 

Laguerre 

A  real,   A  >  0 

Legendre 

Ak  real,   Ak  =^A,A>0 

Jacobi 

Ak  real,   Ak  -  k  +  1 

Kautz 

Ak  complex,   Ak=AKk+jAlk,     ARk>0 
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3.4  GLFF  Kernel  Derivations  and  Properties 

The  Laguerre,  Legendre,  Jacobi,  and  Kautz  GLFF  networks  have  the  interesting 
property  that  their  kernels  are  orthonormal.  Hence,  we  call  these  orthonormal 
generalized  linear  feedforward  (OGLFF)  networks.  Next  these  kernels  will  be 
derived  from  first  principles.  This  is  done  in  the  continuous-time  domain.  The 
derivations  will  make  evident  their  relationship  to  causal  LTI  system.  The  connection 
between  the  network  structures  and  the  mathematical  description  of  causal  LTI 
systems  is  the  primary  reason  they  are  able  to  approximate  these  systems  so  well. 

3.4.1  Continuous-Time  Laguerre  Kernel 

Consider  the  sequence  U/U)"^  \.    Regarding  this  as  the  solution  to  an  ordinary 

differential  equation  that  represents  a  causal  LTI  system,  the  system  would  consist  of 
one  pole  with  multiplicity  n+1.  Orthogonalizing  this  sequence  yields  a  new  sequence 
{/„(/, A)}  given  by 


f„\ 


ln(t,A)  =  ^Ue-A'Y 


(-2Atj 


<=o  V  l  J 


i\ 


,/2>0,/L>0,r>0  (3.11) 


Now  taking  the  Laplace  transform  gives  the  Laguerre  model  kernel 

Ln(s,A)  =  I  {Zn(U)}  =  ^f-^41  (3.12) 

s-  A  \s  +  A  J 


3.4.2  Continuous-Time  Kautz  Kernel  (Complex  Time  Scale  Vector) 

Consider  the  sequence  {ex"'}  where  the  A„'s  are  distinct  complex  scalars  (i.e., 
components  of  the  time  scale  vector  A).   Since  the  poles  are  complex  they  must  exist 


31 


in  conjugate  form  in  order  to  satisfy  the  readability  requirement.    Orthogonalizing 
this  sequence  yields  a  new  sequence  {kn(t, A )}  given  by 

^U)  =  JTRe{A,„^''},      n>0,  t  >  0,  Re{/i,}  >  0,  V/  (3.13) 

It  is  now  shown  how  to  obtain  the  coefficients  AIU.   Taking  the  Laplace  transform  of 
the  new  sequence  gives  the  Kautz  model  kernel 


K2  iU^l^ja^"^^^^^^'^  n-l 


^,,,+i(^)  =  M^,,,+l(^)r- .  i    w  =  0,l ■ 


(3.14) 


Express  these  functions  in  partial  fraction  form.    Grouping  each  pair  of  time  scale 
elements  (At  ,A\ )  gives  the  kernel  expression 


K„(U)  =  2 


2  Aii      ,  .  2  Aw 


S  +  A.        S  +  X. 


,     n  -  2m  and  n  -  2m  +  1 


(3.15) 


ii. 


Hence  the  coefficients  /\,„  are  twice  the  residue  at  the  pole  -A,  for  the  ;/    tap  kernel 


3.4.3  Continuous-Time  Kautz  Kernel  (Real  Time  Scale  Vector) 

Consider  the  sequence   [cjA"']   where   the  i„'s   are  distinct   real  scalars   (i.e., 

components  of  the  time  scale  vector  A).    Orthogonalizing  this  sequence  results  in  a 
new  sequence  {ki:(t,A  )}  given  by 


kn(t,A)  =  £  A,/  ^  n>  0,  t  >  (U,.  >  0,  Vi 


(3.16) 


i=0 


It  is  now  shown  how  to  obtain  the  coefficients  Ani.   Taking  the  Laplace  transform  of 
the  new  sequence  gives  the  Kautz  model  kernel 
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^(a)=/R(a)}=^n£^L=S^j-  (3-17) 

s  +  An%ts  +  Ai      l=Qs  +  Al 
where  the  last  expression  is  simply  the  partial  fraction  expansion  of  the  product  form 
of  Kn(t,  A).  Hence  Ani  is  the  residue  of  the  pole  at  -A\  for  the  n    tap  kernel.   It  will  be 
real  since  the  time  domain  function  kn(t,  A)  is  real. 

3.4.4  Continuous-Time  Legendre  and  Jacobi  Kernels 

Consider  again  the  sequence  [e*"'}  where  the  /l„'s  are  distinct  real  scalars.  This  is 

the  same  sequence  from  which  the  Kautz  kernel  with  a  real  time  scale  vector  is 
derived.  This  kernel  is  described  by  equation  (3.17).  Now  consider  the  kernel  with  a 

"special"  time  scale  parameter.  Namely,  if  A  -  (AVA-,,...,AN)  with  Ak  =2^LA,  A>0 
then  we  obtain  the  Legendre  kernel.  Likewise,  if  Ak  =k  +  l  we  obtain  the  Jacobi 
kernel. 

3.5  Impulse  and  Step  Response  Modeling 

The  models  discussed  in  this  chapter  can  accurately  represent  any  signal  or  system 
belonging  to  L  (K).  We  will  be  most  interested  in  modeling  causal  LTI  systems  using 
sampled  step  response  data.  In  this  case,  a  step  function  is  applied  to  a  continuous- 
time  system  to  generate  its  step  response.  The  step  response  is  sampled  to  obtain  a 
discrete-time  representation  of  the  system.  Optimal  continuous-time  and  discrete- 
time  GLFF  models  are  determined  using  this  data  (optimization  will  be  discussed  in 
Chapter  4).  For  example  using  an  N+l  tap  Kautz  filter  the  approximating  transfer 
function  is 
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ff(z)  =  H(z)  =  5>A(z,A)  (3.18) 


n=0 


Suppose  the  discrete-time  impulse  response  is  desired.   How  can  it  be  computed? 
Three  simple  methods  are  now  given. 

1.  Using  the  discrete-time  model   H(z),  pass  the  impulse  x(k)=d(k)  through  this 
model  to  generate  an  approximation  to  h(k). 

2.  Since  x(t)  is  a  unit  step  then  X(s)  =  l/s.   If  g(t)  is  the  step  response  of  the  system 
then  g(t)  =  h{t)  *  x(t)  or  G(.v)  =  H{s)  I  s.  So 

H  ( s)  =  sG(  s)  =  D(  s)G(s)  (3.19) 

or 

H(z)=D(s)G(s)\     .,  =  D(z)G(z)  where  D(c)  =  -—  (3.20) 

'-ttrr  r  z + 1 

Using  the  sampled  step  response  data  g(k),  pass  this  through  D(c)  to  generate  h(k). 

3.  Since  h(t)  =  — —  use  a  numerical  approximation  to  the  derivative  (backward 

dt 

difference  perhaps)  to  obtain  h(k).  In  this  case 


dRit) 


dt 


„g(kT)-g((k-\)T)  =  g(k)-g(k-\):im 
t=a  '  ' 


3.6  Conclusion  (Model  Formulation) 

This  chapter  has  dealt  with  structures  of  models  useful  for  approximating  causal 
LTI  systems.  The  generic  model  M(p)  was  stated.  It  was  shown  how  this  formulation 
can  be  used  to  represent  FIR,  IIR,  and  GLFF  networks.  Several  GLFF  kernels  were 
derived,  and  discussion  was  provided  to  illustrate  their  connection  to  the  assumptions 
of  the  systems  under  study.  Tables  of  discrete-time  and  continuous-time  GLFF  model 
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kernels  were  included.  These  kernels  are  the  most  popular  in  the  literature  and  as  the 
derivations  of  the  Laguerre,  Kautz,  Legendre,  and  Jacobi  revealed  they  are  not 
theoretical  abstractions,  but  orthogonalized  versions  of  solutions  to  differential 
equation  representations  of  several  important  classes  of  systems.  Now  that  the  model 
structure  M(p)  has  been  defined,  optimization  of  the  parameter  set,  p,  becomes  the 
focus. 


CHAPTER  4 
OPTIMAL  PARAMETER  SELECTION 


In  this  chapter  the  selection  of  the  optimal  (or  near  optimal)  parameter  set  for 
GLFF  networks  is  discussed.  Recall  p  =  {N,  w,A  }  where  TV  is  the  dimension  vector, 
w  is  the  weight  vector,  and  A  is  the  time  scale  vector.  Optimization  for  these  three 
vectors  is  handled  separately  below.  In  fact,  the  ability  to  separate  the  parameters  and 
their  optimization  techniques  is  another  benefit  of  the  GLFF  formulation. 

4.1  Model  Dimension 
Consider  a  GLFF  model  consisting  of  a  set  of  time  domain  basis  functions 
{<t>„(t,A)}'n=Q  with  frequency  representation  {<t>n(s,A )}'  0-    The  model  structure  is  a 
weighted  sum  of  these  basis  functions  (kernels).  In  the  s-domain 

M(p)  =  ^p-  =  £  u'„0„(.v,/l)  =  vv70(.a  )  (4.1) 


where 


w  =  (w0,w1,...,wiV)T 

O(.s',/l)-(O0(^,A),cD1(.v,/l),...,Ov(.v,A))r 


(4.2) 


In  the  time  domain 


v(0 
where 


5>A(U) 


71=0 


*x(t)  =  [wT0(t,A)]*x(t)  (4.3) 
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</)(t,A)  =  (</)0(t,Al<p](t,A),...,<pN(t,A))T  (4.4) 

In  this  thesis,  the  theory  will  focus  on  modeling  the  impulse  response  of  a 

continuous-time  system.     Given  sampled  step  response  data,  continuous-time  and 

discrete-time  GLFF  networks  will  be  constructed  to  approximate  the  desired  system. 

Hence  the  input,  x(t),  will  be  considered  an  impulse  and  the  desired  signal,  d(t),  will 

be  the  system  impulse  response. 

The  selection  of  N  is  now  discussed.    N  is  a  scalar  for  the  GLFF  model  class. 

Suppose  N  is  a  fixed  finite  positive  integer.  Then  d(t)  can  be  approximated  by 

d(0=y(0  =  |>,AM)  (4.5) 

4.1.1  Completeness 

To  ensure  that  y(t)  — >  d(t)  as  N  — >  °°  the  property  of  completeness  is  employed.  A 
set  of  functions  {0„(^)}J=O  is  called  complete  [Su71]  on  the  space  L2([a,b])  if  there 

exists  no  function  d(t)eL2([a,b])  such  that 

b 

jd(t)(/)n(t,A)dt  =  0,       V/i  =  0,1,2,...  (4.6) 

a 

An  equivalent  definition  is  as  follows.     Given  a  piecewise  continuous  function 
d(t)e  L2([a,b])  and  an  £  >  0  ,  there  exists  an  integer  N>0  such  that 

b 

\(d(t)-y(t)fdt<£  (4.7) 
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where  y(t)  is  given  by  the  equation  (4.5).  All  of  the  GLFF  kernel  sets  are  complete  on 

[0,°o).     An  example  of  an  incomplete  set  of  orthonormal  functions  {0„(^)}J=O 

defined  on  [0,27r)  is 

0„(r,/l)  =  0„(O  =  sin(m)  (4.8) 

The   function   d(t)=cos(mt)   is   nonzero   on    [0,2/r)    however   it    is   orthogonal    to 

|^(/,/l)}J=oon  this  interval.  Namely, 

In 

$cos(mt)0n(t,A)dt  =  O,  Vn,m  (4.9) 

0 

4.1.2  Orthogonal  Projection 

When  [<pn(t,A  )}     is  an  orthonormal  set  then 

v(0  =  |>,A(U)  (4.10) 


n=0 


spans  an  N+l  -dimensional  Hilbert  Space,  S  (see  Appendix  B).  The  Orthogonal 
Projection  Theorem  states  that  for  every  d{t)  in  a  Hilbert  space  H,  it  can  be 
decomposed  as  d(t)  =  y(t)  +  e(t)  where  v(/)eS,  e(t)eS  .  S  is  the  orthogonal 
complement  of  S  and  H  =  S  ©  S  .  So  y(t)  can  be  considered  as  the  orthogonal 
projection  of  d(t)E  H  onto  the  space  S.  Namely 

v(0  =  ProJs(</(0)  (4.11) 

and  {<pri(t,A)}      is  an  orthonormal  basis  of  S. 
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4.1.3  "Optimal"  Selection  of  N 

Completeness  is  important  because  it  guarantees  that  if  the  approximation  of  d(t) 

using  N+l  weighted  terms  of  {</>n(t,A)}      is  not  within  a  prescribed  tolerance  then 

more  terms  can  be  included  to  eventually  reach  the  tolerance  requirement.  Increasing 
N  may  not  improve  the  approximation  of  d(t)  when  a  set  of  incomplete  functions  is 
used. 

Having  the  completeness  property  provides  the  following  choices  when  a  GLFF 
network  realization  for  a  given  TV  is  inadequate: 

1.  Increase  N  until  the  prescribed  tolerance  is  reached. 

2.  Select  other  kernels  from  the  GLFF  set  and  examine  the  approximation  error. 

For  a  given  N,  if  none  of  the  GLFF  kernels  chosen  provide  satisfactory  results 
then  N  must  be  increased.  If  this  is  done  using  several  kernels  the  rate  of  convergence 
for  each  kernel  as  N  is  increased  could  be  monitored.  Choose  from  this  set  the  kernel 
that  yields  the  minimum  N  needed  to  achieve  the  required  tolerance.  Other  options 
are  to  either  raise  the  tolerance  level  or  choose  a  different  performance  criteria.  One 
suggestion  for  setting  the  tolerance  is  to  choose  the  acceptable  error  power  to  be  a 
certain  percent  (5%  perhaps)  of  the  system  power. 

Davila  and  Chiang  [Dav95]  suggest  another  approach  to  estimating  the  dimension 
vector  N  =  (N\,  A^)7  for  a  general  IIR  structure.  This  method  examines  the 
eigenvectors  of  the  input/output  data  covariance  matrix.  When  the  model  orders  are 
overestimated,  zeros  will  appear  in  the  noise  subspace  eigenvectors.  The  number  of 
zeros  found  can  be  used  to  estimate  the  model  orders  (dimension  vector  components). 
This  procedure  could  be  applied  to  GLFF  networks. 
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4.2  Weight  Vector 


One  of  the  greatest  advantages  to  working  with  a  feedforward  model  is  that  the 
calculation  of  the  optimal  weight  vector  w  (in  the  squared  error  sense)  is  simple.  This 
is  true  even  for  GLFF  structures.  Although  the  realizability  requirement  for  GLFF 
networks  demands  the  weight  vector  be  real,  the  optimal  solution  for  w  is  now 

derived  assuming  w,  d(t),  and  {</>n(t,A)}      are  generally  complex.  Realizable  optimal 

solutions  are  obtained  by  setting  the  imaginary  components  equal  to  zero. 

4.2.1  Optimal  Selection  of  vv 

Assume  cl{t)  is  the  impulse  response  of  the  desired  system  (i.e.,  d{t)  =  hit))  and 
approximate  it  with  y(t)  given  by 

,v 

>'(0  =  5>,A(U)  (4.12) 

/i=0 

The  parameter  set  p  =  {N,w,A}  is  to  be  chosen  so  the  error  e(t)  =  d{t)  -  y(t)  is 
minimum  in  the  integrated  squared  error  sense.  Namely,  choose  p  to  minimize 

J  =  j\d(t)-y(tfdt 

0 

oo  ft  oo 

=  j\d(tfdt  -  £  w]d\t)^{uX  )dt  (4. 13) 

o  '=o       o 

-|)  w; ]dit)f,it,A  )dt  +  X|wJ2 ]\0,(t,A  )fdt 

1=0  o  ,=0  o 

If  vv„  =  a„  +  jfr„  then  set =  0  and =  0  to  find  vv„  that  minimizes  J.  Hence 

dan  dbn 

-\      j  oo  oo  oo 

—  =  -J  d\t)<t>n{t,X)dt  -jd(t)fn(t,A  )dt  +  2anj\0n(t,A,)\2dt  =  0       (4.14) 

'<  0  0  0 
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so 


JRe[d(t)fn(t,A)]dt 

an=^~z (4.15) 

j\</>n(t,Afdt 

0 

When  {0n(?,A)}      is  an  orthonormal  set,  the  bottom  integral  equals  unity  yielding 

an=\Re[d(t)fn(t,A)]dt  (4.16) 

o 

Likewise  bn  can  be  derived 
dJ 


dK 


j\d(t)fn(t,A)dt-j\dt(t)(/))i(t,A)dt  +  2bn^n(t,Afdt^0      (4.17) 


so 


]lm[d(t)fn(t,A)]dt 
K=*—z (4-18) 

)\(/>n(t,A)\2dt 


Again,  if  {0n(t,A)}      is  orthonormal  then  equation  (4.18)  reduces  to 

bn=jlm[d(t)fn(t,A)]dt  (4.19) 

0 

Hence,  combining  the  above  results  gives 

wn=an  +  jbn=jd(t)fn(tMdt  (4.20) 

0 

When  {</>n(;t,A)\     is  real  then  w„  is  real  and  can  be  written  as 

wn=]d(t)^n(t,X)dt  (4.21) 
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So 


\v  =  jd(t)0(t,A)dt 


(4.22) 


is  the  weight  vector  that  minimizes  the  energy  of  e(t).   From  this  point  forward  only 
real  weight  vectors  will  be  considered. 

4.2.2  Relationship  to  Wiener-Hopf  Equations 

It  is  easy  to  show  that  w  is  in  fact  the  solution  to  the  Wiener-Hopf  (or  Normal) 
equations.  The  performance  function  can  be  expanded  as  follows 

J  =  \e\i)dt  =  \{d{f)-  y(t)fdt  =\(d(t)-wT0(t,A ))" dt 


=  fd2(t)dt-2wTjd(t)0n(t,A)dt  +  wT  j<f>(t,A)<f)T(t,A)dt 


\v 


now  define  the  following  terms 


de  =  \d\t)dt 

0 

p  =  \d(t)</)(t,A)dt 

0 

R  =  J0(t,A)f(t,A)dt 


(4.23) 


(4.24) 


Then  the  performance  function  J  becomes 

J  -  d  -  2wT p  +  w  Rw 


(4.25) 


To  minimize  this  function  take  the  gradient  of  J  with  respect  to  the  weight  vector  w 
and  equate  to  zero 


Vh,7  =-2p  +  2Rw  =  0 


(4.26) 


or 
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Rw  =  p  (4.27) 

These  are  known  as  the  Wiener-Hopf  or  Normal  equations.  Solving  for  w  gives  the 
optimal  solution 

w  =  R~1p  (4.28) 

Now  the  OGLFF  kernels  {</>n(t,A)}     form  an  orthonormal  set  so 

R  =  j<t>(t,A)f(t,A)dt  =  I  (4.29) 

0 

Hence  the  weight  vector  becomes 

w  =  R~xp  =  p=  (d(t)0(t,A)dt  (4.30) 

o 

which  is  the  solution  derived  in  the  previous  section. 

4.2.3  Minimum  Squared  Error  Performance  Function 

Now  that  the  optimal  weight  vector  has  been  derived,  the  minimum  squared  error 
can  be  computed.  Recall 

J  =  de-2wTp  +  wTRw  (4.31) 

For  the  optimal  weights  it  was  found  that  (using  OGLFF  kernels)  w  =  p  and  R  =  I. 
Plugging  these  into  the  equation  for  J  yields 

Jopt  =de-  wTw  (4.32) 

The  above  optimal  weight  vector  w  and  minimum  squared  error  Jopt  was  derived  using 
a  fixed  time  scale  A.  Notice  that  w  and  Jopt  are  functions  of  A.  Hence  Jop!  can  be 
written 

Jopl(A)  =  de-wT(A)w(A)  (4.33) 
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to  illustrate  the  importance  of  A  on  the  optimal  solution.     Appendix  C  presents  a 
derivation  of  equation  (4.32)  in  the  ^-domain. 

4.2.4  Computing  w  Using  Step  Response  Data 

In  Chapter  5  applications  of  the  GLFF  networks  are  given.  The  systems  there  are 
continuous-time  systems.  The  available  information  of  these  systems  is  sampled  step 
response  data.  GLFF  networks  will  be  used  to  approximate  the  input/output  step 
response  behavior.  Methods  for  computing  the  optimal  weight  vector  w  using 
sampled  step  response  data  are  now  given. 

Least  Squares  Solution.  [Hay91]  Consider  the  model  dimension  N  and  time  scale 
vector  A  fixed  or  previously  chosen  to  be  optimal.  Discussion  of  model  dimension 
was  given  in  Chapter  4.1.  Discussion  of  the  optimal  time  scale  will  be  given  in 
Chapter  4.3. 

Since  the  information  is  sampled  (discrete-time)  step  response  data  the  following 
definitions  are  needed.  Let  the  input,  x{t),  be  the  unit  step  and  g(t)  be  the  response  of 
the  system  under  question.  If  the  step  response  is  sampled  at  interval  T  then  define 
g(k)  =  g(kT),  &=0,1,  ...,  AT  as  the  desired  signal.  Define  also 

xi(k)  =  klh  sample  at  the  /,h  tap  of  the  GLFF  network  (4.34) 

y(k)  =  klh  output  of  the  GLFF  network 
In  this  discrete-time  case,  the  matrix  R  and  vector/?  are  given  by 

R  =  [r]  whcrtrv=(xi,xj)  =  ^xi(k)xj(k)  (4.35) 

and 
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KPn) 


where     pi  =  {g,xi)  =  ^g(k)xi(k)  (4.36) 


*=o 


This  gives  the  normal  equations  Rw  =  p  with  solution  w  =  R'lp.  This  is  the  Wiener- 
Hopf  solution.  Notice  in  this  case  that  R  will  not  be  the  identity  matrix.  Only  when 
the  input  is  an  impulse  (or  white  noise)  will  R  =  I.  Nevertheless,  the  Wiener-Hopf 
solution  provides  the  optimal  weight  vector  (for  fixed  TV  and  X)  regardless  of  the 
input. 

Integral  Equation  Solution.  [Clu91]  A  more  direct  evaluation  of  the  weight  vector 
solution  given  in  Chapter  4.2.1  is  now  considered.  Recall 

wn=\d(t)<pn(t,A)dt,     n  =  0,\,...,N  (4.37) 

o 

is  the  optimal  weight  solution  assuming  an  orthonormal  basis.     For  the  system 

modeling  problem  this  is  accomplished  by  applying  an  impulse  (or  white  noise)  into 

the  system.  In  which  case  d(t)=h(t),  the  impulse  response.  Consider,  for  example,  the 

Laguerre  kernels  given  by  [0n(t,Z)}  =0={^(^^)}  =0-  Some  of  its  properties  can  be 
utilized  to  derive  a  simplified  integral  equation  solution.  Again  work  with  g(k),  k=0,l 
,...,  K  samples  of  the  continuous-time  step  response  g{t).  Since  the  continuous-time 
impulse  response  is  h(t)  =  dg(t)/dt,  the  optimal  weights  can  be  written  as 


wn  =  ]h(t)Ut,A)dt 


dgit) 


i^-ln(t,A)dt  (4.38) 


o 


dt 


=(s«)/„(a)|;-j*(o^U 
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where  integration  by  parts  has  been  used.    For  a  stable  system  g(t)  it  must  have 
|g(oo)|  <  oo.  Also  using  the  property  of  the  Laguerre  functions  [Par71] 


-*  dn 


Ut,A  )  =  Vll— ^—  (t ne'1M)  (4.39) 

yields 

lim/„(r,A)  =  0  (4.40) 

Hence 

<?/„(U) 


„=-U/)^4^  (4.41) 

i  dt 


0 


Let  r„  be  the  steady  state  response  time  of  g(t);  namely,  Tss  is  the  time  which  for  all 
t  >  Tss  the  output  can  be  consider  constant,  g(t)  =  gss  =  constant.  In  this  case  the 
integral  can  be  broken  up  as 

dln(U) A     .  }dln(t,A) 


r         01  u, A  r0iu,.. 

0  T„ 

«  dt 


°  (4.42) 

dln(t,A) 


From  the  Laguerre  function  property  above  it  is  easy  to  compute  the  derivative 
relationship 

^4^  =  -A/,,(/,A)-2A£/,(U)  (4.43) 

dt 


1=0 


The  weight  computation  becomes 

n-\  T, 

SJ 

i=0    0 


wn  =  2A^\g(t)li(t,A )dt  +  A\g(t)ln(t,A  )dt  +  gJn(Ts,A)  (4.44) 
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Up  to  now  the  continuous-time  step  response  g(t)  has  been  assumed.  To  compute  w„ 
using  the  sampled  step  response  g(k)  an  integral  approximation  technique  must  be 
used  such  as  Euler  integration 

]g{t)ll{t,A)dt^^g{k)ll(k,A)  (4.45) 


t=o 


where 

T 

A  =  — —     (integration  interval)  (4.46) 

K 

Finally  this  gives  the  integral  equation  weight  solution 

wn=2AA^j^g(k)li(k,A)+AA^g(k)ln(k,A)  +  gJn(Tss,A)  (4.47) 


i"=0  *=0  t=0 


4.2.5  Adaptation  of  w 

Note  that  the  weights  can  be  computed  by  solving  a  set  of  linear  (Wiener-Hopf) 
equations.  This  linear  relationship  among  the  weights  corresponds  to  a  convex 
performance  function  with  one  minimum  point.  So  the  weight  vector  can  be  easily 
solved  in  an  adaptive  manner  using  a  gradient  search  technique  such  as  Newton's 
method  or  LMS  [Wid85].  Unfortunately,  this  nice  relationship  does  not  hold  for  the 
time  scale  vector.  It  is  (in  general)  non-convex  with  multiple  stationary  points.  The 
number  of  which  depends  upon  the  dimension,  N,  of  the  model  used. 

4.3  Time  Scale  Vector 

Optimization  of  the  parameter  set  elements  N  and  w  have  been  discussed  to  this 
point.    Their  solutions  are  fairly  easily  found  for  GLFF  networks  and  are  the  only 
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elements  of  p  that  are  required  for  FIR  and  HR  models.  As  has  been  previously 
stated,  the  power  of  GLFF  networks  is  in  local  feedback.  This  feedback  is  controlled 
by  the  time  scale  vector.  Although  IIR  structures  also  have  feedback,  it  can  not  be 
separated  and  optimized  independently  from  the  tap-to-tap  weights  like  the  GLFF 
structures. 

Unfortunately,  the  elements  of  the  time  scale  vector  are  not  linearly  related;  hence, 
the  performance  function  J  is  usually  non-convex  w.r.t.  these  parameters.  As  such, 
the  bulk  of  the  work  in  finding  an  optimal  p  involves  finding  an  optimal  A.  This 
section  is  devoted  to  this  issue.  It  begins  by  examining  several  methods  that  currently 
exist  for  either  determining  exactly  or  approximately  the  optimal  A.  Each  method 
usually  addresses  a  particular  GLFF  kernel,  most  often  the  Laguerre.  A  few  papers 
have  considered  the  non-orthogonal  Gamma  and  a  few  discuss  optimization  of  the 
more  complex  Kautz  kernel.  Although  the  exact  optimal  solution  is  sought,  Kautz 
has  pointed  out  that  non-optimal  solutions  (which  may  be  derived  in  a  simplified 
manner)  can  still  produce  excellent  results.  His  comments  suggest  that  one  should 
not  search  to  exhaustion  for  optimality  if  it  is  possible  to  get  close  with  much  less 
work. 

After  the  overview  of  optimization  methods  found  in  the  literature,  we  propose  a 
new  technique  relying  heavily  on  complex  variable  theory.  This  is  a  method  that 
extends  the  concepts  of  squared  error  analysis  to  the  complex  domain.  Under  certain 
assumptions,  the  optimal  time  scale  can  be  fairly  easily  determined  by  simply  looking 
for  zero  crossings. 
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4.3.1  Optimal  Selection  of  /I  (Historical  Review) 

In  this  section  a  summary  is  given  of  current  methods  found  in  the  literature  to 
either  locate  or  approximate  the  optimal  time  scale  A.  Each  method  usually  addresses 
a  particular  GLFF  kernel,  most  often  the  Laguerre.  A  few  papers  have  considered  the 
non-orthogonal  Gamma  and  a  few  discuss  optimization  of  the  more  complex  Kautz 
kernel. 

The  methods  presented  in  the  literature  can  be  grouped  into  the  following 
categories  according  to  the  employed  strategy  or  desired  objective: 

1.  Approximation  via  the  impulse  response. 

2.  Asymptotic  solutions. 

3.  Achieving  prescribed  measurement  objectives  (other  than  minimum  squared 
error). 

4.  Matched  filter  approach. 

5.  Moment  matching. 

6.  Satisfying  necessary  conditions  for  stationarity  of  the  squared  error  function. 

7.  Iterative  and  adaptive  search  techniques. 

Approximation  via  the  impulse  response.  These  methods  approximate  A  given  a 
priori  knowledge  about  the  impulse  response  h(t)  of  the  system.  The  required 
information  is  different  for  each  approach.  Either  h(t)  is  needed  in  analytic  form  or 
tabular  form  (samples  of  h(t)  at  discrete  points,  t  =  tk,  k  =  0,1,2,  ...  ).  These  methods 
appeal  to  the  necessary  requirement  for  all  stable  and  realizable  linear  systems; 
namely,  its  transfer  function  H(s)  must  be  representable  as  a  rational  function  with  no 
right-half  plane  poles  or  multiple  j-axis  poles.  Under  this  assumption,  h(t)  can  be 
represented  as  a  sum  of  exponentials  of  the  type 

(A0  +  A,t+...+Aj'")eAt'  (4.48) 

where 
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**=*#  + Am    ^^0  (4-49) 

A  frequency  domain  method  for  approximating  A  is  as  follows  [Su  71].  Suppose 
h{t)  is  known  (or  approximation  to  it)  and  its  Laplace  transform  H{s)  is  determined. 
Note  that  the  form  of  h{t)  that  is  given  may  not  produce  a  rational  H(s).  Now  find  a 
rational  approximation  Ha(s).  For  instance  use  the  expression  for  H(s)  and  expand  it 
into  a  Taylor  series  about  s  =  0  to  obtain 

H(s)  =  a0+ais  +  a2s2  +•••  (4.50) 

Then  a  rational  function  can  be  found  to  approximate  H(s).  One  way  to  accomplish 
this  is  to  use  a  Pade  approximant.  The  poles  of  this  approximation  can  be  used  as  an 
initial  guess  for  the  poles  of  H(s).  These  pole  locations  could  be  incrementally 
adjusted  to  locally  optimize. 

A  time  domain  approach  is  to  use  point  matching.  The  basis  of  this  technique  is 
attributed  to  Prony  [Ham73].  The  strategy  is  to  force  an  approximate  impulse 
response  h(t)  to  match  h(t)  at  evenly  spaced  points.  Assume  h(t)  is  known  at  points 
separated  by  the  spacing  At.  Using  the  approximation 

1=1 
it  is  possible  to  create  N  linear  equations  (using  the  known  points  of  h(t))  whose 
solution  gives  coefficients  a0,  fli.  •••>  tf/v-i  of  an  A^h  order  characteristic  polynomial 

xN  +  aN_]xN~l  +  aN_2xN"2  +  — h  axx  +  a0  =  0  (4.52) 

Once  the  coefficients  are  known,  find  the  roots  of  this  polynomial  x\,  X2,  ...,  xN.  They 
are  related  to  the  time  scale  parameters  by 
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At=—\nx,  (4.53) 

'     At 

Another  time  domain  approach  involves  the  use  of  an  orthogonal  signal  set 
[Kau54].  Assume  h(t)  and  its  first  N  derivatives  are  known.  Using  these  functions, 
generate  a  set  of  orthogonal  functions  {0„(O}„=O  satisfying  the  relationships 

</>Q(t)  =  h(t) 

^(t)  =  h(t)  +  bw(/)0(t) 

<p2(t)  =  ii(t)  +  b2]</)l(t)  +  b20(/)0(t)  (4.54) 

0N(t)=hiN)(t)  +  bNiN_l)</>N_l(t)  +  ---+bN](/)](t)  +  bNO</)o(t) 

where  the  &,„„'s  can  be  found  such  that 

]</>j(t)<f>k(t)dt  =  Sjk  (4.55) 

0 

th  th 

Since  each  $,(0  is  a  linear  combination  of  the  derivatives  (from  the  0  to  the  n  )  of 
h(t)  then  $,(r)  can  be  rewritten  in  terms  of  Iv'if),  j  =  0,1 ,...,/?— 1 .  Replace  each  h^it) 
by  A!  to  get  the  characteristic  equation 

#  +  ^-l)^  +  ^(n-2)^"2  +'  "+  M  +  «„0  =  0  (4-56) 

The  roots  of  these  equations  give  the  time  scale  parameters  Ani,  i=l,2,...,n  for  model 
order  n  for  each  n  =  1,2,. ..,N. 

Asymptotic  solutions.  Several  methods  are  offered  that  are  based  on  power  series 
representations  of  the  Laguerre  model.  Schetzen  [Sch70],  [Sch71]  presents  an 
analytic  method  for  determining  the  asymptotic  optimal  time  scale  /L  where 

/L  =  lim  AN  (4.57) 

N— »~ 
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This  is  based  on  the  region  of  convergence  (ROC)  of  the  power  series  equivalence  of 
the  Laguerre  model.  Justification  for  /L,  in  the  place  of  An  is  made  on  the  basis  that 
most  engineering  applications  require  the  use  of  a  large  number  of  terms  TV  in  the 
Laguerre  expansion  to  achieve  accurate  results.  In  this  case,  /L  is  viewed  as  a  good 
approximation  to  the  optimal  time  scale  An. 

Bruni  [Bru64]  offers  another  approach  involving  optimization  around  the  ROC  of 
a  Laguerre  model  expansion.  Assuming  h(t)  is  the  impulse  response  of  an  LTI 
lumped  parameter  system 

M 

ft&)  =  £V*      pl  =  of  +  j«iffll>0  (4.58) 

i=0 

decompose  a  Laguerre  model  representation  of  h(t)  into  M  infinite  series 

oo  oo         M 

MO  =  £H;/,,(U)=X5>m/,,(U)  (4.59) 


n=0  n=0  1=1 


•th 


where  vvm  is  the  n  coefficient  of  the  expansion  of  the  i  exponential  in  h(f).  With 
this  formulation  the  Fourier  transform  of  h{t)  can  be  represented  as  a  sum  of  M 
infinite  complex  geometric  series.  M,  gives  the  magnitude  of  the  term  in  the  i 
infinite  series.  M,  should  be  minimized  for  each  i  in  order  to  achieve  the  most  rapid 
convergence  of  the  series.  Since  M,  is  a  function  of  A,  there  is  an  optimal  A  (denoted 
Ai)  for  each  M,.  Choose  as  the  optimal  A  for  the  entire  model,  the  A{  that  minimizes 
the  largest  M,.  Methods  similar  to  this  one  are  also  given  in  [War54]  and  [Hea55]. 

Achieving  prescribed  measurement  objectives.  In  addition  to  locating  A  that 
minimizes  the  squared  error,  approaches  that  minimize  other  measurement  criteria 
have  also  been  considered.    One  such  method  employed  when  using  the  Laguerre 
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signal  set  {i„(t,A)}n=Qis  given  in  [Par71].  Consider  the  approximation  of  d{t)  using 
an  N+ 1  tap  Laguerre  model 

d(t)=y(t)  =  fjwJn(t,A)  (4.60) 

Now  define  the  following  quantities 

oo  oo 

\td2(t)dt                         \td2{t)dt 
Mx=^ and     M2  =  -^ (4.61) 

jd2(t)dt  jd\t)dt 

o  o 

M\  is  a  measure  of  how  rapidly  d(t)  decays.  Mi  is  a  measure  of  how  smoothly  d(t) 
decays.     Parks  show  that  the  time  scale  A  given  by  A  =  JMjM2  minimizes  the 

maximum  error  over  the  class  of  all  signals  with  these  measurements.  The  advantage 
to  this  approach  is  that  a  complete  knowledge  of  the  signal  to  be  approximated  is  not 
required. 

In  [Fu93]  a  similar  approach  is  taken  when  using  the  discrete-time  Laguerre 
model.  In  this  case,  the  optimal  A  is  found  that  minimizes  the  measurement 

./  =  f>w;  (4.62) 

n=0 

This  performance  measure  is  chosen  because  it  linearly  increases  the  cost  of  each 
additional  term  used  in  the  Laguerre  expansion.  This  enforces  a  fast  convergence 
rate.  The  optimal  A  is  obtained  using  a  discrete-time  formulation  of  M\  and  Mi 
above.  A  relationship  for  J  is  derived  that  depends  only  on  M\,  Mi,  and  A.  The 
optimal  A  is  then  determined  to  be 
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A  = ■    '  (4.63) 

2M,  - 1  +  ^4MlM2  -  Ml  -2M2 

Matched  filter  approach.  In  [Kuo94]  a  matched  filter  bank  is  used  to  estimate  the 

optimal  A  for  the  Gamma  model.  Recall  that  {gn(t,A )}      is  a  non-orthogonal  GLFF 

signal  set.  Let  d(t)  be  approximated  by 

rf(0sy(0  =  5>.ft.(U)  (4-64) 

A  matched  filter  is  constructed  based  on  the  optimal  A  condition  (g^+l(t,A),d(t))  =  0 

where  gf,+l(t,A)  is  the  component  of  gN+\U,  A)  that  is  orthogonal  to  the  previous  basis 
components.  A  bank  of  such  filters  is  created  where  each  bank  uses  a  different  A.  To 
optimize  A  over  the  region  (0,  Amax)  an  M-bank  filter  structure  is  created  where  the  mth 
filter  uses  the  time  scale  A  =  Amax  I  M.  When  the  desired  signal  d(t)  is  applied  to  the 
filter  bank,  a  change  in  sign  of  the  M  outputs  will  occur  between  filter  m  and  in+l  for 
a  value  of  A  that  satisfies  the  optimality  condition.  When  there  are  multiple 
occurrences,  the  associated  A  for  each  case  must  be  used  to  determine  which  one 
achieves  the  minimum  squared  error. 

Moment  matching.    The  method  presented  in  [Cel95]  uses  a  Gamma  model  to 
approximate  a  signal  d{t).  The  Gamma  model  is  conceptualized  as  a  rotating  manifold 

in  space.  The  kernels  {g,M^)}n=>  are  tne  basis  vectors  that  define  the  N+\ 
dimensional  space  of  operation  and  the  time  scale  parameter  A  is  the  feature  that 
allows  for  rotation.  The  optimal  A  is  estimated  using  a  moment  matching  method. 
This  method  attempts  to  estimate  A  by  equating  a  finite  number  of  time  moments  of 


54 

y(t)  (the  output  of  the  Gamma  model)  to  those  of  d(t).  Namely,  it  is  required  that 
mn(y(t))  =  mn(d(t)),    n  =  0,...,N  where 

mn(y{t))  =  \tny{t)dt  (4.65) 

0 

Using  these  functions  an  N+\  order  polynomial  in  A,  p(A),  is  formulated.  The  roots  of 
p(A)  are  estimates  of  the  local  minima  of  the  squared  error  performance  function  J. 
The  root  that  minimizes  J  is  chosen  as  the  optimal  A. 

Satisfying  necessary  conditions  for  stationarity  of  the  squared  error.  Some  of  the 
most  popular  methods  for  optimizing  the  time  scale  of  the  Laguerre  model  center 
around  the  necessary  conditions  for  stationarity  of  the  squared  error  performance 
function  J.  For  instance,  approximating  d{t)  with  an  N+\  tap  Laguerre  network 

d{t)^wn{A)ln(t,A)  (4.66) 


n=0 


where  the  weight  vector  components  wn(A)  have  been  given  an  argument  in  A  to 
emphasis  the  dependence  of  these  weights  on  the  time  scale.  For  a  fixed  A,  the 
squared  error  is  minimized  when  the  weights  are  computed  as 

wn(A)  =  ]d(t)ln(t,A)dt  (4.67) 

0 

The  squared  error  function  reduces  to 

J  =  ]d\t)dt-%w2n(A)  (4.68) 

0  "=0 

To  minimize  J  (over  A)  the  evaluation  of  the  X  must  be  maximized.  This  requires  that 

^Evv,;(i)  =  X2w,?a)-^vv„a)  =  0  (4.69) 

dA  „=0  „=o  dA 
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which  yields  the  interesting  relation 

wN(A)wN+l(A)  =  0  (4.70) 

A  proof  of  this  assertion  is  given  in  Appendix  D.  Hence  the  optimal  value  of  A  is 
either  a  root  of  wN(A)  =  0  or  wN+i(A)  =  0.  Solving  for  A  analytically  is  a  tedious  task 
because  w^(A)  =  0  is  a  polynomial  in  A  of  degree  (N+M)  if  d(t)  is  a  system  with  M 
poles.  The  first  proof  of  this  (using  mathematical  induction)  was  given  by  Clowes 
[Clo65]  and  was  restricted  to  the  case  where  d{t)  had  only  simple  poles.  A  more 
concise  and  less  restrictive  proof  was  given  in  [Kin69].  This  paper  also  derives 
expressions  for  the  2nd  derivative 
'N  +  h 


,2     N 


^— [(# +l)wfl+lU)-NwN+lU)wN_1(AJ\         forwN(A)  =  0 

2A  (4.71) 

^[{N  +  2)Wn+2(A)wn(A)-(N  +  l)wl(A)}    forwN+1(A)  =  0 

2/t 

These  equations  may  be  used  to  determine  whether  a  stationary  point  is  a  maximum 
or  minimum.  In  [Mas90]  a  discrete-time  formulation  of  the  stationary  conditions  is 
formulated  and  the  authors  use  numerical  techniques  to  find  the  roots  of  the 
polynomials  wN(A)  =  0  and  wN+\(A)  =  0.  [Sil93]  derives  these  conditions  for  the 
Gamma  model.  The  derivation  given  by  Clowes  was  valid  only  for  systems  with 
rational  transfer  functions.  [Wan94b]  extends  the  results  to  approximating  any  L"(R) 
stable  linear  system.      [Sil95]      describes  another  efficient  way  to  compute  the 

derivatives  — -wn(A)  and  uses  this  to  compute  high-order  approximations  of  the 
dA 

weights  wN(A)  and  wN+J(A)  via  Taylor  series  and  Pade  approximants.  This  procedure 

is  illustrated  in  both  continuous-time  and  discrete-time  when  approximating  a  system 

excited  by  an  arbitrary  input. 
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Iterative  and  adaptive  search  techniques.  When  representing  a  system  with  a 
discrete-time  GLFF  network,  the  elements  of  the  time  scale  vector  are  restricted  to 
bounded  regions  to  guarantee  stability.  In  particular,  the  discrete-time  Laguerre 
model  requires  Ae  (—1,1).  One  method  of  optimization  is  to  employ  exhaustive  search 
over  this  interval,  seeking  the  A  that  minimizes  the  squared  error  performance 
function.  Since  this  would  be  done  in  discrete  steps,  the  only  way  to  get  continually 
closer  to  the  optimal  A  is  to  increase  the  number  of  points  in  the  search  region. 
[Kin77]  and  [Mai85]  discuss  minimizing  the  squared  error  by  employing  a  Fibonacci 
search  technique  over  A.  This  is  an  interval  elimination  procedure  wherein  the  region 
in  which  the  minimum  lies  is  sequentially  reduced  using  a  predetermined  number  of 
steps. 

[Nur87]  uses  an  adaptive  technique  whereby  the  optimal  A  is  chosen  to  be  the 
center  of  gravity  of  the  presumed  poles  of  the  system.  The  method  works  as  follows. 
Initialize  A,  say  Aq.  Using  the  Laguerre  model,  with  this  A,  the  parameters  of  a 
difference  equation  representation  of  the  system  are  estimated.  From  the 
characteristic  equation  determine  the  poles  of  the  estimated  system.  Choose  as  Aj, 
{=1,2,...  the  center  of  gravity  of  these  poles.    Repeat  this  process  until  \Al+l  -  At\<  S 

where  <5"is  a  small  number.  Once  the  center  of  gravity  converges  to  a  fixed  point  stop 
the  adaptation. 

4.3.2  Optimal  Selection  of  A  (Complex  Error  Analysis) 

A  new  approach  to  estimate  the  optimal  time  scale  is  now  offered.  The 
methodology  used  here  is  one  of  extending  the  squared  error  performance  function  to 
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the  complex  domain.  This  is  accomplished  by  defining  a  new  complex  performance 
function.  It  will  be  shown  how  this  function  is  derived,  how  it  relates  to  the 
conventional  squared  error  function,  and  what  additional  information  it  offers.  The 
optimal  time  scale  can  be  determined  by  looking  for  zero  crossings  of  the  imaginary 
component  of  the  new  complex  performance  function. 

For  systems  consisting  of  a  cascade  of  2nd  order  subsystems  containing  pairs  of 
complex  conjugate  poles,  examples  will  be  presented  that  illustrate  how  the  imaginary 
component  of  the  new  complex  performance  function  can  be  used  to  obtain  the 
optimal  (in  the  minimum  squared  error  sense)  time  scale  parameter.  These  examples 
will  solve  the  system  identification  problem  with  a  GLFF  model  containing  the 
Laguerre  kernel. 

4.3.2.1  Derivation  of  the  complex  performance  function.  In  Appendix  C  an  s- 
domain  description  of  the  squared  error  performance  function  is  discussed.  It  is  given 
by 

J  =  ^E(s)E(-s)ds  (4.72) 

c 

where  C  is  a  contour  taken  in  the  CCW  direction  enclosing  the  entire  left-hand  plane 
(LHP)  in  s,  E(s)  =  D(s)-  wTQ>(s,X),  and  w  is  a  real  weight  vector.  For  a  fixed  time 
scale  X  the  optimal  weight  vector  is 

w  =  ^iD(s)®(-s,A)ds  (4.73) 

c 

yielding  the  performance  function 

j  =  de  -wTw  (4.74) 

where 
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de=^JD(s)D(-s)ds  (4.75) 

c 

Now  consider  restricting  evaluation  of  the  performance  function  to  the  upper  LHP 

only.  Namely,  create  the  complex  performance  function 

Ju^^JE(s)E(-s)ds  (4.76) 

Make  the  following  definitions: 

due  =  2^  <j>  D{s)D{-s)ds  (complex  energy) 

c"  (4.77) 

wu  =  j-  <J> D(s)<f>(-s,A)ds  (complex  weight  vector) 

c„ 

The  contour  Cu  is  taken  in  the  CCW  direction  enclosing  the  upper  LHP  in  s.  Now  if 
F(s)  is  a  function  which  has  real  and/or  complex  conjugate  pair  poles  then 


^JF(s)ds  =  2Re 


^§F(s)ds 


(4.78) 


hence 

w  =  2Re[w,  1 

r    \  (4.79) 

<*e  =  2ReK] 

By  making  the  above  definitions,  the  problem  of  locating  the  optimal  A  using  Ju  is 
studied.  By  operating  in  the  complex  space  there  is  information  that  can  be  used  to 
offer  additional  insight  into  the  optimal  solution  of  the  time  scale  parameter  A.  This 
information  is  not  available  using  squared  error  analysis  because  the  imaginary 
component  of  Ju  is  canceled  when  calculating  J.  This  is  due  to  the  contribution  from 
the  poles  of  D(s)  in  the  lower  LHP  in  s.  The  squared  error  performance  function  can 
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be  considered  as  a  mapping  onto  R  and  the  complex  performance  function  as  a 
mapping  onto  C. 

A  simplified  expression  for  /„  is  derived  as  follows 

Ju=^JE(s)E(-s)ds 

c, 

=  -^  §[D(s)-  wt$>(s,A)][d(-s)  -  wT®(-s,A)]ds 

=  ^§D(s)D(-s)ds-wT^-JD(s)<l>(-s,A)ds 
cu  c, 

- wT -^§D(-s)Q>(s,A  )ds  +  wT ■£$ §®(s,A  )$>T(-s,A)dsw 


(4.80) 


Using  the  definition  of  the  complex  energy,  due,   and  complex  weight  vector,  wu,  as 
well  as  the  following  relations 


vv7  -^&(&(s,A)<$>T(-s,A)dsw  =  ±wTw 


(4.81) 


and 


^JD(-s)0(s,A,)ds  =  iw 


(4.82) 


gives 


Note  that 


-  d„-  wTw„ 


ue  it 


2Re[7u]  =  2Re[rfHJ  -  v/2Re[wJ  =  de  -  wTw  =  J 
which  is  what  is  expected  since 


J  =^§E(s)E(-s)ds  =  2Re 


-^JE(s)E(-s)ds 


(4.83) 


(4.84) 


=  2Re[7H]  (4.85) 
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4.3.2.2  Complex  performance  function  examples  and  experimental  results.  Three 
test  cases  have  been  chosen  to  study  the  complex  performance  function,  Ju.  One  2nd 
order  system  and  two  4  order  systems  are  modeled  using  a  GLFF  network  with  a 
Laguerre  kernel.  The  transfer  function  for  each  example  follows: 


Example  1:  (2nd  order  system,  rapidly  decaying  impulse  response) 


H(s)  = + (4.86) 

s+2-j     s+2+j 


Example  2:  (4(  order  system,  rapidly  decaying  impulse  response  with  slight 
oscillation) 


H(s)  = ! + ! + ! + ! (4.87) 

s  +  2-j     s  +  2  +  j     s  +  l-2j     s  +  l  +  2j 


Example  3:  (4n  order  system,  slowly  decaying  impulse  response  with  oscillation) 


H(s)  = + + - + 1 (4.88) 

5  +  0.4-;     s  +  0.4  + j     5  +  0.6-2.8;     s  +  0.6  +  2.8; 


Experimental  results  for  each  example  are  tabulated  and  graphed  on  the  next  few 
pages.  In  Figures  4-1  and  4-2  the  squared  error  performance  function,  J  (solid  curve), 
and  twice  the  imaginary  component  of  the  complex  performance  function,  2Jm[Ju] 
(dash  curve),  are  plotted  for  the  first  four  model  orders  for  example  1,  Figures  4-3 
and  4-4  for  example  2,  and  Figures  4-5  and  4-6  for  example  3.  Notice  the  imaginary 
component   has   zero  crossings   near  the   stationary  points   of  the   squared  error 
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performance  function.  Im[7„]  crosses  negative-to-positive  near  local  minima  of  J.  So 
the  complex  performance  function  provides  an  estimate  of  the  optimal  time  scale 
through  the  imaginary  component.  In  fact  Im[/J  can  be  considered  an  estimate  of 
V  XJ(A)  that  is  best  near  Aopt. 

Numerical  results  are  provided  in  Tables  4-1,  4-2,  and  4-3  to  illustrate  how  close 
the  nearest  zero  crossing  (nearest  ZCU)  of  2Im[7„]  is  to  the  optimal  time  scale  as  the 
model  order  increases.  Also  the  squared  error  is  computed  at  the  optimal,  J(Aop!),  and 
nearest  ZCU,  J(Anzcu),  time  scales.  For  this  analysis,  Im[/J  was  scaled  by  2  to  provide 
similar  scaling  to  J  since  J=2Re[Ju]\  however,  the  zero  crossings  are  unaffected  by 
this  scaling.  By  examining  J(Aopt)  and  J(AnZcu)  it  is  found  that  when  using  a  Laguerre 
network  with  sufficient  order  to  approximate  the  unknown  system,  no  degradation 
occurs  by  choosing  A  at  the  nearest  ZCU.  In  fact  the  small  difference  in  Aop,  and  AnZcu 
is  partly  due  to  numerical  approximation  error  when  searching  for  minima  of  J  and 
roots  of  Im [/,<]. 
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2  3 

Time  Scale,  X 

Figure  4-1:7  for  example  1  (orders  =  0  to  3). 
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Figure  4-2:  J  and  2Im|7,/]  for  example  1. 
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2  3  4 

Time  Scale,  X 

Figure  4-3:  J  for  example  2  (orders  =  0  to  7). 
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Figure  4-4:  J  and  2Im[7w]  for  example  2. 
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2  3  4 

Time  Scale,  A. 


Figure  4-5:  J  for  example  3  (orders  0  to  11). 
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Figure  4-6:  J  and  2Im[7(/]  for  example  3. 
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Table  4-1:  Performance  function  characteristics  for  example  1 


N 

A„p, 

AiZCu 

J(Aopt)  in  dB 

JUnzcu)  in  dB 

0 

2.47 

2.48 

-24.7 

-24.6 

1 

1.60 

1.58 

-33.8 

-33.7 

2 

2.31 

2.31 

-49.9 

-49.9 

3 

1.92 

1.92 

-60.7 

-60.7 

Table  4-2:  Performance  function  characteristics  for  example  2. 


N 

Aft 

A„ZCu 

J(A,opt)  in  dB 

J(A„ZOl)  in  dB 

0 

2.99 

3.12 

-11.6 

-11.6 

1 

1.50 

1.43 

-15.6 

-15.5 

2 

2.58 

2.60 

-20.8 

-20.8 

3 

1.81 

1.78 

-24.8 

-24.7 

4 

2.45 

2.47 

-29.3 

-29.3 

5 

1.93 

1.92 

-33.3 

-33.3 

6 

2.39 

2.40 

-37.7 

-37.7 

7 

2.00 

2.00 

-41.8 

-41.8 

Table  4-3:  Performance  function  characteristics  for  example  3. 


N 

* 

AiZCu 

J(Aol)<)  in  dB 

J(AnZcu)  in  dB 

0 

0.82 

0.34 

-0.5 

-0.2 

1 

2.30 

2.46 

-2.7 

-2.7 

2 

1.29 

1.21 

-5.2 

-5.1 

3 

2.12 

2.12 

-9.9 

-9.9 

4 

2.87 

2.89 

-12.7 

-12.7 

5 

3.56 

3.59 

-14.1 

-14.1 

6 

4.21 

4.22 

-14.7 

-14.7 

7 

4.82 

4.76 

-14.9 

-14.9 

8 

2.24 

2.24 

-17.1 

-17.1 

9 

2.62 

2.64 

-19.6 

-19.5 

10 

2.25 

2.24 

-22.9 

-22.9 

11 

2.57 

2.56 

-25.5 

-25.5 
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4.3.3  Further  Analysis  of  the  Complex  Performance  Function 

In  the  previous  section,  it  was  shown  that  the  squared  error  and  complex 
performance  functions  can  be  expressed  as 


7". 


J  =  d—w  w 


h   =  Clue  -  WT™U 


(4.89) 


where 


de=i-§D(s)D(-s)ds       dtte=^$D(s)D(-s)ds 
w  =  ^JD(s)®(-s,A)ds    wu  =  ^§D(sm-s,A)ds 


(4.90) 


and 


D(s)  =  desired  system 
<J>(s,/l)  =  L(s,A)  the  Laguerre  kernel 

C  =  contour  enclosing  the  entire  left  -  hand  plane  (LHP)  in  s 
Cu  =  contour  enclosing  the  upper  left  -  hand  plane  (LHP)  in  s 


It  was  also  shown  that 


^  =  2ReK] 
w  =  2Re[wJ 


(4.91) 


Writing  J  and  Ju  in  normalized  form  (relative  to  the  total  energy  of  the  system)  gives 

(4.92) 


J=(de-wTw)/de 
lt-K-wTwu)/de 


Decomposing  due  and  wu  into  their  complex  component  forms 

due  =  dueR  +  Jduel 
W„   =  WuR  +  JWyl 


(4.93) 


yielding 


de  =  2Ridue]  =  2duel 
w  =  2Re[wJ  =  2w„R 


(4.94) 
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So  we  can  write  J  and  Ju  in  terms  of  components  of  the  complex  energy  signal  and 
complex  weight  vector. 

J  =  (2dut,R-(2wJ(2WuR))/(2dlieR) 


"  (4.96) 

=  (12due-wTuRwu)/du 


lueR 


Now  we  are  interested  in  2Im[/u],  twice  the  imaginary  component  of  the  complex 
performance  function.  This  has  the  expression 

2lm[ji]  =  (duel-2wTuRwul)/dueR  (4.97) 

This  extension  to  the  complex  domain  has  shed  light  on  new  information.  The 
imaginary  component  of  Ju  has  been  shown  to  have  roots  (or  zero  crossings)  near  the 
stationary  points  of  J,  one  of  which  provides  an  approximation  to  the  optimal  time 
scale  solution. 

In  an  effort  to  capitalize  on  this  discovery,  the  analytic  performance  function,  Ja, 
will  be  defined  and  studied  as  a  method  of  achieving  the  same  information  while 
requiring  less  a  priori  knowledge  of  the  unknown  system  D(s).  As  will  be  seen,  Ja  is 
composed  of  an  energy  signal  and  a  weight  vector  that  is  the  analytic  signal 
representation  of  the  associated  squared  error  components. 

4.3.3.1  Analytic  performance  function.  Consider  the  components  of  J„.  The 
complex  energy  signal,  duc,  can  be  written  as 
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due=^JD(s)D(-s)ds 

0  °° 

--i-\D{a)D{-(j)d(j  +  i-\D{ja))D{-ja))dco 

-^o  0 

~  r      o 

^JDUo))D(-ja))dQ)  -j  £  J  D(cy)D{-cj)d(7 

o  J      L    ■— 

Likewise,  the  complex  weight  vector,  wu,  can  be  written  as 
wu=-^$D(sm-s,A)ds 

0  oo 

=  i-\D((j)<S>(-(j,A)d(j  +  ^\D(jco)<$>(-ja),A)da) 

-oo  0 

Now  collecting  real  and  imaginary  components  gives 

~  r      o 

■^JD(JQ))^(-jO),X)do)  -Re^JD(am-a,A)dcr 

o  J        L    -~ 

i     r    ° 

^JD(ja))®(-ja),A)dct)  -Im  £  j  D((T)®(-cr,A)d<T 


(4.98) 


(4.99) 


w.  =  < 


Re 


>  + 


Im 


(4.100) 


Consider  the  components  of  wu  along  the  +jco  axis  only.    This  is  equivalent  to  the 
spectrum  of  the  analytic  signal  constructed  from  w.    Namely,  for  a  signal  /(/),  the 


analytic  signal  is  specified  by 


f.{t)*f(t)  +  jf(t) 


(4.101) 


where    f(t)  -  H\f(t)\,   H  is  the  Hilbert  Transform.     The   spectrum   of   fa{t)    is 
equivalent  to  the  spectrum  of  /(?)  for  0  <  co  <  °°  and  zero  for  -°°  <  co  <  0.    Now  we 


define  the  following  functions: 
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dae  =  [real  part  of  dM  along  +  jco  axis]  +  j [imaginary  part  of  due  along  +  jco  axis] 

+  j[0]  (4.102) 


-±JD(jco)D(-jco)dco 

o 

±JD(jco)D(-jco)dco 


and 


wfl  =  [real  part  of  wu  along  +  jco  axis]  +  7 [imaginary  part  of  wu  along  +  jco  axis] 


Re 


^jDO^)O(-yftU)^ 

0 

---^JD(jcom-jco,A)dco 


+  j\lm 


-^^D(jco)<P(-jco,A)dco 


(4.103) 


Jfle  is  called  the  analytic  energy  signal  and  wa  is  called  the  analytic  weight  vector.  As 
in  the  case  of  the  complex  performance  function  we  define  the  analytic  performance 


function,  Ja, 


or  in  normalized  form 


J„=d„  -  w  vv 


a  ae 


Ja={dae-wTwa)/de 

Now  decomposing  the  components  of  Ja  into  their  complex  form  gives 

dae  ~  daeR  +  Jdae, 
Wa   =  W„R  +  JWa, 


(4.104) 


(4.105) 


(4.106) 


where 


laeR 


-^JD(jco)D(-jco)dco 


(4.107) 


d«,  =  0 


and 
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WaR  =  Re 


wal  =  Im 


^JD(j(o)^(-jo},A)dQ) 

0 

£JD(ja))®(-jct),A,)da) 


(4.108) 


The  goal  is  to  rewrite  J  and  Ja  in  terms  of  analytic  energy  signal  and  analytic  weight 
vector  components  only.  Recall 

w  =  i-JD(s)^(-s,A)ds 

c 

=  £  J D(jco)®(-jco,A)dco 

o  ~ 

=  i  JD(jco)^(-jco,A)dct)  +  ^JD(ja})0(-jQ),A)dQ) 

-oo  0 

-^JD(jco)<3>(-jco,A)d(o    +  ^\D(jco)^(-jco,A)dco  (4.109) 

0 

=  2Re  ^JD(jco)0(-jco,A)dct) 

0 

=  2Re[wJ 


=  2w 


aR 


Likewise 


J.  =  2Re 


-^JD(jco)D(-jco)dco 


=  2Re[<*J 


=  2d 


aefl 


So  we  can  write  J  and  7a  in  terms  of  components  of  dae  and  wa  as  follows 

y  =  (2^-(2wrt/?)r(2wfl/?))/(2^) 

Ja={dae-{2waR)Twa)/(2daeR) 

=  {2dae-wTaRWa)/daeR 


(4.110) 


(4.111) 


(4.112) 
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We  are  interested  in   2Im[/J,  twice  the  imaginary  component  of  the  analytic 


performance  function.  This  has  expression 

2Jm[ja] 


-2w[Dw 


K 


aR     al  I      aeR 


(4.113) 


4.3.3.2  Analytic  performance  function  examples  and  experimental  results.  Again 
consider  the  three  examples  used  to  study  the  complex  performance  function.  Figures 
4-7,  4-8,  and  4-9  compare  the  complex  and  analytic  performance  functions.  For  each 
of  the  three  examples  the  squared  error  performance  function,  J  (solid  curve),  twice 
the  imaginary  component  of  the  complex  performance  function,  2Im[7,J  (dash  curve), 
and  twice  the  imaginary  component  of  the  analytic  performance  function,  2Im[7a] 
(dot-dash  curve),  are  plotted  for  the  first  four  model  orders.  Numerical  results  are 
provided  in  Tables  4-5,  4-6,  and  4-7  to  illustrate  how  close  the  nearest  zero  crossing 
(nearest  ZCU)  of  2Im[7u]  and  the  nearest  zero  crossing  (nearest  ZCa)  of  2Im[/fl]  are  to 
the  optimal  time  scale  as  the  model  order  increases.  Also  the  squared  error  is 
computed  at  the  optimal,  J(Aopl),  nearest  ZCU,  J(A„zcu),  and  nearest  ZCa,  J(A,nzca),  time 
scales. 

It  is  interesting  to  compare  Ju  and  Ja  for  these  three  examples.  The  difference  in 
these  functions  is  that  Ju  includes  the  line  integral  along  the  negative  craxis.  Consider 
the  complex  energy  due  for  each  example.  From  Table  4-4  we  find  that  duei  is  22%  of 


Table  4-4:  Complex  energy  signal  characteristics. 


Example 

due—  «uefl+  dlie[ 

%  complex 

1 

0.45+0.10/ 

22 

2 

1.98+0.63/ 

32 

3 

9.91+0.31/ 

3 
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dueR  for  example  1,  32%  for  example  2,  and  only  3%  for  example  3.  Observe  from 
Figure  4-9  that  JU(A)  ~  Ja(A),  for  all  X,  for  example  3.  Figures  4-7  and  4-8  show  that 
Ju  and  Ja  are  significantly  different  (but  have  close  zero  crossings  near  A()pl)  for 
examples  1  and  2.  So  the  line  integral  along  the  negative  <j  axis  contributes 
significantly  to  the  evaluation  of  J  for  examples  1  and  2  but  not  example  3.  This 
results  in  the  similarity  of  Ju  and  Ja  for  the  last  case.  Clearly,  the  relative  values  of 
dueR  and  duei  are  application  dependent. 

We  have  shown  that  even  though  Ju  and  Ja  are  different  functions,  their  imaginary 
components  have  roots  that  are  similiar  near  Aopt.  Since  Ja  can  be  computed  directly 
from  a  sampling  of  d{t),  it  could  be  used  to  provide  the  optimal  time  scale  estimate  in 
lieu  of  Ju. 
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Figure  4-7:  J,  2lm[Ju],  and  2Im[7n]  for  example  1. 


Table  4-5:  Performance  function  characteristics  for  example  1. 


N 

A()pt 

AnZCu 

A„ZCa 

J(Alip,)  in  dB 

J(/\nZCll)  in  dB 

JU„ZCa)  in  dB 

0 

2.47 

2.48 

2.48 

-24.7 

-24.6 

-24.6 

1 

1.60 

1.58 

1.49 

-33.8 

-33.7 

-31.1 

2 

2.31 

2.31 

2.29 

-49.9 

-49.9 

-49.5 

3 

1.92 

1.92 

2.62 

-60.7 

-60.7 

-54.6 
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Figure  4-8:  J,  2Im[7,J,  and  2Im[/fl]  for  example  2. 


Table  4-6:  Performance  function  characteristics  for  example  2. 


N 

Ajp/ 

AnZCu 

AiZOi 

J(Al>pf)  in  dB 

JUn/n,)  indB 

J(A„ZCn)  in  dB 

0 

2.99 

3.12 

3.08 

-11.6 

-11.6 

-11.6 

1 

1.50 

1.43 

1.31 

-15.6 

-15.5 

-14.8 

2 

2.58 

2.60 

2.19 

-20.8 

-20.8 

-19.0 

3 

1.81 

1.78 

3.16 

-24.8 

-24.7 

-22.5 

4 

2.45 

2.47 

2.40 

-29.3 

-29.3 

-29.2 

5 

1.93 

1.92 

1.86 

-33.3 

-33.3 

-32.8 

6 

2.39 

2.40 

2.14 

-37.7 

-37.7 

-34.8 

7 

2.00 

2.00 

1.96 

-41.8 

-41.8 

-41.5 
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Figure  4-9:  J,  2Im[7(/],  and  2Im[/fl]  for  example  3. 


Table  4-7:  Performance  function  characteristics  for  example  3. 


N 

■^op/ 

AiZCu 

A>ZGi 

J(A(lfH)  in  dB 

JUnzc)  in  dB 

JU„ZCa)  in  dB 

0 

0.82 

0.34 

0.55 

-0.5 

-0.2 

-0.4 

1 

2.30 

2.46 

2.43 

-2.7 

-2.7 

-2.7 

2 

1.29 

1.21 

1.21 

-5.2 

-5.1 

-5.0 

3 

2.12 

2.12 

2.13 

-9.9 

-9.9 

-9.9 

4 

2.87 

2.89 

2.89 

-12.7 

-12.7 

-12.7 

5 

3.56 

3.59 

3.60 

-14.1 

-14.1 

-14.1 

6 

4.21 

4.22 

4.22 

-14.7 

-14.7 

-14.7 

7 

4.82 

4.76 

4.79 

-14.9 

-14.9 

-14.9 

8 

2.24 

2.24 

2.22 

-17.1 

-17.1 

-17.0 

9 

2.62 

2.64 

2.60 

-19.6 

-19.5 

-19.6 

10 

2.25 

2.24 

2.28 

-22.9 

-22.9 

-22.8 

11 

2.57 

2.56 

2.62 

-25.5 

-25.5 

-25.3 
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4.3.4  Time  Domain  Synthesis  Using  GLFF  Networks 

For  the  discrete-time  setting,  one  of  the  fundamental  limitations  of  an  FIR  model 
is  its  lack  of  memory  depth;  namely,  its  inability  to  approximate  (with  low  order)  a 
system  with  an  infinitely  decaying  impulse  response.  GLFF  networks  do  not  have 
this  problem  due  to  the  presence  of  local  feedback  in  their  structure. 

How  well  a  Laguerre  network  approximates  the  impulse  responses  of  the  three 
example  systems  given  by  equations  (4.86),  (4.87),  and  (4.88)  is  now  illustrated. 
Here  the  system  impulse  response,  d(t),  (h(t)  =  d{i)  since  the  input  is  an  impulse),  the 
Laguerre  network  output  when  using  the  optimal  time  scale  yo(t),  and  the  Laguerre 
network  output  when  using  the  nearest  ZCU  time  scale  yz(t)  are  plotted.  In  example  1 
the  impulse  response  decays  rapidly  and  has  no  significant  oscillation.  yo(t)  and  yz(t) 
are  plotted  in  Figure  4-10  using  a  1st  order  Laguerre  model.  In  example  2  the  impulse 
response  is  also  rapidly  decaying  but  with  a  slight  oscillation.  Plots  of  yo(t)  and  yz(t) 
are  given  in  Figure  4-11  for  a  5th  order  model.  Example  3  is  a  system  having  a  slowly 
decaying  impulse  response  with  a  strong  oscillation.  Plots  of  yo(t)  and  yz(t)  are  given 
in  Figure  4-12  for  a  10th  order  model.  This  oscillation  is  difficult  to  approximate  with 
a  Laguerre  model  since  only  a  single  real  pole  (time  scale)  is  used  in  the 
approximation.  Increasing  the  order  of  the  network  could  further  reduce  the  error. 
Another  option  is  to  use  a  Kautz  model  with  a  pair  of  complex  conjugate  poles. 

For  all  three  examples  we  observe  that  yo(t)  and  yz(t)  are  inseparable.  The  point 
to  be  made  in  this  section  is  that  when  using  A  obtained  from  the  appropriate  zero 


77 


1.5 


0.5 


Example  1:  d(t),yo(t),yz(t) 


12  3  4  5 

Figure  4-10:  Impulse  response  and  approximations  using  a  1st  order  Laguerre  model. 


Example  2:  d(t),yo(t),yz(t) 


Figure  4-11:  Impulse  response  and  approximations  using  a  5l  order  Laguerre  model. 


Example  3:  d(t),yo(t),yz(t) 


Figure  4-12:  Impulse  response  and  approximations  using  a  10l  order  Laguerre  model. 
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crossing  of  Im[7,J  (or  Im[7a])  as  the  time  scale  for  the  model,  the  same  result  can  be 
achieved  as  when  using  A  obtained  by  searching  J  for  the  global  minimum.  Also  the 
squared  error  performance  function  can  always  be  computed  using  J  =  2Re[/J. 

4.3.5  Conclusion  of  the  Complex  Performance  Function 

For  systems  composed  of  sets  of  real  and/or  complex  conjugate  pair  poles,  a  new 
complex  performance  function,  Ju,  was  derived  when  using  a  GLFF  network  as  a 
system  model.  The  imaginary  component  of  Ju  has  zero  crossings  near  the  stationary 
points  of  the  squared  error  performance  function  /.  Experimentally  it  is  determined 
that  as  the  order  of  the  model  is  increased,  the  zero  crossings  of  Im[/„]  near  minima 
of  J  converge  to  these  stationary  points,  one  of  which  will  be  the  optimal  time  scale, 
Aop[.  Rather  than  searching  the  non-convex  function  J  for  a  global  minimum,  locate 
the  roots  of  lm[Ju(A)],  Ar,  and  evaluate  2Re[/„(^r)]-  One  of  these  will  be  close  to  A()p,; 
namely, 

Kt  =  argf  min{2Re[/war)]}j  (4.114) 

where 

Ar  =  { roots  of  lm[Ju(k)] } . 
In  addition,  we  have  defined  the  new  analytic  performance  function  and  derived  its 
closed  form  solution.  It  was  shown  that  J  could  be  written  in  terms  of  components  of 
Ju  and  Ja.  Expressions  for  the  imaginary  components  of  J„  and  Ja  were  also  given. 
The  value  of  these  functions  in  determining  the  optimal  time  scale  of  GLFF  kernels 
was  illustrated  using  three  example  systems.  The  important  equations  from  the  above 
sections  are  now  summarized: 
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Squared  Error  Performance  Function 


J  =  \ 


(de  -  wTw)/de  {Squared  error  components} 

(dlieR  -  2wTuRwuR)/dlieR     {Complex  error  components}  (4. 1 15) 

(daeR  ~  2wTaRwaR)ldaeR     {Analytic  error  components} 


Complex  Performance  Function 

J  =Ud    -wTRw  )/d  R  (4.116) 

u         \2      ue  uR      u  J I       ueR  v  ' 


Analytic  Performance  Function 

Ja={idae-wTaRwa)/daeR  (4.117) 

where 

de=^-§D(s)D(-s)ds  w  =  ^-JD(s)<S>(-s,A)ds 

c  c 

due=i-§D{s)D{-s)ds  wu  =  ^$D(sm-s,A)ds  (4.118) 

dae=ii\D{jco)D{-j(o)d(0  wa=^JDUci))^(-jO),A)dQ) 

0  0 

If  D(s)  is  known  analytically  or  in  closed  form  we  can  employ  the  residue  theorem 
to  compute  de,  due,  dae,  w,  wu,  wa.  If  these  functions  must  be  determined  numerically 
then  the  analytic  performance  function  can  be  easily  computed.  The  only  requirement 
is  a  sampled  data  set  from  an  impulse  or  step  response  test.  In  addition,  J  can  always 
be  computed  from  Ju  or  Ja.  This  illustrates  that  the  squared  error  solution  is 
completely  embedded  within  the  complex  and  analytic  performance  functions. 
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4.4  Conclusion  (Optimal  Parameter  Selection) 
This  chapter  has  addressed  optimization  techniques  for  GLFF  models,  M(p). 
Optimization  of  each  element  of  the  parameter  set  p  -{N,w,A}  is  handled 
separately.  Model  dimension  element  optimization  generally  involves  increasing  TV 
until  a  prescribed  error  tolerance  is  achieved.  Monitoring  the  rate  of  convergence  for 
several  GLFF  kernels  could  also  aid  in  choosing  the  model  yielding  the  smallest  N. 
The  optimal  weight  vector  minimizes  the  integrated  squared  error  and  is  found  by 
solving  a  set  of  Wiener-Hopf  equations.  A  direct  integral  evaluation  method  was  also 
given  that  takes  advantage  of  properties  of  the  model  kernels.  The  weight  vector 
solution  is  a  function  of  the  time  scale  vector  A.  Optimizing  A  is  the  most  difficult 
task  because  the  minimum  squared  error  solution  is  nonlinear  w.r.t.  these  parameters. 
A  number  of  optimization  approaches  exist  in  the  literature  and  they  were  discussed. 
Finally,  a  new  method  was  proposed  that  extends  the  squared  error  analysis  to  the 
complex  domain.  Simply  locating  zero  crossings  and  choosing  the  one  that 
minimizes  J  gives  the  optimal  A.  Examples  were  provided  to  demonstrate  the  theory. 


CHAPTER  5 
APPLICATIONS 


Applications  of  the  theory  presented  in  the  previous  chapters  will  now  be 
discussed.  In  this  thesis,  GLFF  networks  will  be  used  to  approximate  sampled-data 
systems.  In  particular,  the  step  response  of  a  system  used  in  laboratory  test  and 
evaluation  of  guided  weapon  systems  will  be  modeled.  Other  areas  where  GLFF 
models  have  proven  useful  are  first  surveyed. 

The  GLFF  model  class  encompasses  many  structures  that  have  been 
independently  formulated  and  studied  in  the  past.  This  includes  the  tap-delay  line, 
Gamma,  Laguerre,  Legendre,  Jacobi,  and  Kautz  models  defined  in  both  discrete-time 
and  continuous-time  domains.  The  unification  of  these  structures  is  one  contributory 
component  of  this  thesis;  however,  it  is  worth  reviewing  the  previous  history  of  these 
models  to  gain  insight  into  the  many  applications  that  they  have  been  able  to  address. 

Network  synthesis.  The  earliest  work  dealt  primarily  with  orthogonal  model 
structures  in  an  effort  to  synthesize  networks  and  transient  responses  of  linear  systems 
[Arm57],  [Arm59],  [Kau54],  [Hea55],  [Bru64],  [Clo65],  [Kin77],  [War53].  In 
[Lee33]  and  [Hug56]  synthesis  of  electrical  circuits  was  the  focus. 

System  identification  and  modeling.  Similarly  as  the  tools  became  more 
developed  the  nomenclature  and  approach  shifted  to  system  identification  and  system 
modeling  [Bod94],  [Eko94],  [Kin79],  [Lin93],  [Mas90],  [Mas91],  [Nur87],  [Nin94], 
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[OH87],  [OH94],  [Per91],  [Sil94],  [Sil95],  [Van94],  [Wah91],  [Wah94],  [Moh88]. 
The  basic  principles  are  the  same  as  in  network  synthesis  only  now  the  modeling 
effort  involves  more  than  just  the  network  impulse  response.  A  further  refinement 
includes  modeling  non-causal  systems  [Eko95],  nonlinear  biological  systems 
[Mam93],  probability  density  functions,  [Pab92],  [Tia89].  State  space  formulations 
appeared  [Dum86],  [Gus93],  [Heu90],  [Heu95],  and  modeling  using  step  response 
data  [Wan94b]  and  frequency  response  data  [Clu92],  [Wan95]  have  been  considered. 

Function  approximation.  A  representation  of  an  underlying  system  is  not  always 
desired.  Various  GLFF  models  were  also  useful  in  function  approximation,  signal 
representation,  time  series  analysis  [Bro65],  [Cel95],  [Con94],  [Cle82],  [Hwa82], 
[Kor91],  [Mar69],  [Mai85],  [McD68],  [Ros64],  [Sch71],  [Wah93],  [You62b], 
[Yen62].  Particular  examples  include  representation  of  seismic  data  [Dea64]  and 
electrocardiograms  [You63]. 

Speech.  Applications  of  GLFF  models  used  in  speech  synthesis/analysis  [Man63] 
and  speech  compression  [A1J94]  have  appeared. 

Control  theory.  GLFF  models  are  also  useful  in  control  theory.  Examples  include 
adaptive  controllers  [Dum86],  [Dum90],  [Zer88a],  adaptive  predictor  controllers 
[Els94],  [Guo94],  optimal  control  [Hag91],  [Nur81],  robust  control  [Clu91]  and 
robust  stability  tests  [Aga92]. 

Stochastic  processes.  Stochastic  concepts  have  been  addressed  including 
approximating  cross  correlations  [den93a],  correlations  and  power  density  spectra 
[Lee67],  and  probability  density  functions  [Tia89],  [Pab92]. 
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(Adaptive)  filter  theory.  GLFF  models  have  been  used  as  filters  [Kab94] 
including  adaptive  filters  [den93a],  [den94],  [Per88],  [Pri93].  Systematic  procedures 
to  design  FIR  filters  have  been  introduced  [Fri81],  [Mas90],  [Mas91].  Recursive 
identification  of  time  varying  signals  [Gun90],  2D  filtering  [Tsi94]  and  hardware 
implementations  have  also  been  studied  [Ji95]. 

Other  applications.  Other  examples  where  GLFF  networks  have  been  used 
include  solutions  to  boundary  value  problems  [Ira92],  separation,  enhancement,  and 
tracking  multiple  sinusoids  [Ko94],  and  approximation  of  time  delays  [Lam94], 
[Par91],  [Moh95]. 

Although  the  list  of  applications  above  is  fairly  extensive,  it  is  by  no  means 
exhaustive.  Nevertheless,  the  applicability  of  the  GLFF  modeling  theory  developed  in 
this  research  is  clearly  illustrated. 

5.1  Examples  of  GLFF  Networks  Using  Simulated  Systems 

The  application  focus  of  this  thesis  is  system  modeling;  namely,  it  is  desired  to 
approximate  (model)  LTI  systems  using  GLFF  networks.  The  objective  is  to  obtain 
an  approximation  to  the  system  impulse  response.  Quite  often  this  is  accomplished 
by  optimizing  the  parameters  of  the  model  using  step  response  data.  In  this  section 
several  GLFF  network  realizations  are  constructed  that  approximate  simulated  LTI 
systems.  Optimization  is  performed  using  the  step  response  of  the  desired  systems. 
Real  (non-simulated)  systems  and  their  associated  step  response  data  will  be 
considered  in  Chapters  5.3  and  5.4. 
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Three  simulated  systems  will  be  studied:  1st  order,  2nd  order,  and  4th  order. 

Example  1:  (1st  order) 

#(*)  =  —!-,    r  =  l  (5.1) 

TS  +  i 

Example  2:  (2nd  order) 

H(s)  =  — ^ r,  co„  =  1,£  =  0.4  (underdamped)  (5.2) 

.s  +  2^fyn5  +  ft;; 

Example  3:  (4th  order) 

H(s)=Hi(s)H2(s)  (5.3) 

where 


*W  =    -     „/""' r^iu  =  1>£  =  °-4  (underdamped) 


ffl? 


//,(s)  =  -= — T~^ni  -  2,^2  ~  0-25  (underdamped) 


(5.4) 


In  each  case,  the  step  response  is  generated  and  sampled.  The  sampled  step 
response  data  is  modeled  with  Laguerre,  Gamma,  and  Legendre  models.  The  results 
are  given  below  in  Figures  5-1  through  5-6.  In  each  case,  the  optimal  time  scale 
parameter  was  determined  and  used  to  compute  and  plot  the  step  response  and 
impulse  response  for  6-tap  (N  =  5)  networks.  Figures  5-1,  5-2,  and  5-3  compare  the 
system  and  GLFF  network  step  responses.  Figures  5-4,  5-5,  and  5-6  show  the 
resulting  impulse  responses.  The  optimal  time  scale  parameters  of  the  6-tap  models 
and  associated  minimum  squared  error,  J(Aopt),  are  provided  in  Tables  5-1,  5-2,  5-3. 
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Figure  5-1:  1st  order  step  response  and  model  data. 
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Figure  5-2:  2     order  step  response  and  model  data. 
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Figure  5-3:  4    order  step  response  and  model  data. 
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Figure  5-4:  1st  order  impulse  response  and  model  data. 
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Figure  5-5:  2     order  impulse  response  and  model  data. 
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Figure  5-6:  4th  order  impulse  response  and  model  data. 
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Table  5-1:  1st  order  model  results. 


Model 

N 

A-opt 

J(Aopr)  in  dB 

FIR 

101 

- 

-13.7 

Laguerre 

5 

0.88 

-57.1 

Gamma 

5 

0.05 

-oo 

Legendre 

5 

i          „005< 

^opt,   =  e 

-56.0 

Table  5-2:  2nd  order  model  results. 


Model 

N 

Aopt 

J(Aopt)  in  dB 

FIR 

101 

- 

-13.0 

Laguerre 

5 

0.95 

-57.1 

Gamma 

5 

0.08 

-40.7 

Legendre 

5 

i                     0.03; 
*opt,   =  e 

-37.5 

Table  5-3:  4    order  model  results. 


Model 

N 

A-opi 

J(Aopt)  in  dB 

FIR 

101 

- 

-10.9 

Laguerre 

5 

0.89 

-37.7 

Gamma 

5 

0.09 

-31.2 

Legendre 

5 

i                    0.03/ 

Ku  =  e 

-28.8 
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The  three  examples  considered  above  are  typical  of  the  subsystem  makeup  of 
guided  weapon  systems.  With  fairly  low  order  models  (N=5)  the  GLFF  networks 
have  demonstrated  (in  the  simulated  examples)  they  can  approximate  these  systems 
well.  This  includes  systems  having  slowly  decaying  impulse  responses  and 
oscillatory  modes.  The  most  difficult  task  of  the  modeling  process  is  locating  the 
optimal  time  scale.  This  was  done  using  a  numerical  search  technique  in  the  above 
examples.  For  a  given  model,  if  the  resulting  error  is  too  large  the  model  order  can  be 
increased.  This  process  can  be  continued  until  the  prescribed  error  tolerance  is 
reached. 

5.2  Laboratory  Testing  of  Guided  Weapon  Systems 

Guided  weapon  systems  are  as  important  in  a  military  tactical  mission  as  the 
aircraft  and  ships  that  carry  them.  They  are  the  tools  with  which  a  pilot  can  carry  out 
both  offensive  and  defensive  plans.  As  such,  their  capability  and  reliability  must  be 
well  understood  so  that  they  can  be  operated  within  the  envelope  for  which  they  were 
developed.  To  ensure  their  success  it  is  necessary  that  the  operational  envelope  is 
understood.  This  knowledge  is  acquired  through  extensive  test  and  evaluation  of 
these  systems.  The  objective  is  to  produce  a  smart  dependable  weapon  that  operates 
as  expected  when  launched  from  an  aircraft  or  ship. 

Test  and  evaluation  (T&E)  of  guided  weapon  subsystems  can  be  broken  down 
into  three  phases  [Eic89].  This  breakdown  is  typical  when  a  new  weapon  or 
modification  of  an  existing  weapon  is  under  consideration.  These  phases  are 
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1.  Research/Concept  Exploration 

a)  analytic  study 

b)  open-loop  laboratory  and  field  tests 

2.  Simulation 

a)  software  simulation  and  modeling 

b)  HITL  simulation 

3.  Flight  Test 

a)  captive  carry 

b)  launch 

Research  and  concept  exploration  are  the  very  beginning  of  the  T&E  effort. 
Proposals  must  be  carefully  considered  to  determine  which  (current)  weapons  could 
possibly  meet  the  objective  being  considered.  Can  an  existing  weapon  be  modified  or 
does  it  call  for  a  new  design?  All  these  scenarios  must  be  analytically  explored. 
Flight  tests  are  the  final  hurdle.  It  is  here  where  the  weapon  system's  performance  is 
measured  against  its  expectation.  Before  flight  testing  however,  many  problems  can 
be  worked  out  by  evaluating  a  weapon  system  in  a  simulated  environment  created  in 
the  laboratory.  There  are  many  methods  for  evaluating  the  performance  of 
components  of  a  weapon  system  such  as  parametric  tests  and  open-loop  tests. 
Another  class  of  tests  studies  the  behavior  of  these  systems  in  a  closed-loop 
configuration  and  involves  a  simulation  of  the  real-world  environment  and  the 
weapon's  flight  in  this  environment.  This  class  of  tests  can  be  broken  down  into  two 
types:  all-digital  and  hardware-in-the-loop  (HITL).  In  an  all-digital  closed-loop 
flight,  the  entire  weapon  system  is  simulated  mathematically.  If  the  weapon 
subsystem  models  are  accurate,  all-digital  simulations  can  give  us  a  good  idea  of  how 
the  weapon  will  perform  as  well  as  what  its  limitations  and  vulnerabilities  are.  Once 
subsystems  have  been  hardware  fabricated,  these  can  replace  the  models  in  the 
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closed-loop  simulation.     Ultimately  it  is  desired  to  have  as  much  of  the  weapon 
hardware  in  the  loop  as  possible. 

5.2.1  Objectives 

Several  reasons  why  laboratory  HITL  testing  is  an  integral  part  of  a  weapons 
development  test  program  are  now  given. 

1.  Asset  Availability 

During  the  development  of  a  weapon  system,  often  answers  are  needed  to  key 
questions  about  the  system  to  determine  the  limits  of  its  capability.  Rather 
than  fabricating  an  entire  system  and  conducting  a  flight  test,  many  times 
these  questions  can  be  answered  by  using  simple  prototypes  or  only  subsystem 
components.  This  precludes  having  to  completely  fabricate  a  system  and 
conduct  open  air  flight  tests.  Flight  tests  are  expensive  and  consume  lots  of 
resources  (fully  developed  weapon  system,  aircraft/ships,  open-air  range 
time/space,  money).  Here  is  where  laboratory  testing  offers  a  great  advantage. 

2.  Planning  and  Analysis  Requirements 

As  a  missile  system  reaches  development  stage  and  can  be  flight  tested,  valid 
test  conditions  (test  matrix)  must  be  formulated.  Having  laboratory  results 
(all-digital  and/or  HITL)  as  a  guide,  the  test  matrix  can  be  formulated  with 
intelligence  and  unrealistic  scenarios  can  be  avoided.  Laboratory  tests  can  not 
only  help  define  the  launch  envelope  of  a  weapon,  but  they  can  be  used  to 
predict  its  performance  when  operating  under  valid  conditions. 

3.  Repeatability 

Another  advantage  to  the  laboratory  environment  is  repeatability  testing. 
Once  a  weapon  system  simulation  has  been  developed,  it  can  be  used  to 
conduct  hundreds  and  thousands  of  tests.  Weapon  performance  when 
launched  under  varying  release  conditions  and  "what  if  scenarios  can  be 
examined.  The  laboratory  environment  is  a  major  contributor  to  this  kind  of 
testing  because  real  assets  do  not  have  to  be  expended  to  conduct  many  tests. 
In  "live  fire"  testing  assets  can  usually  be  used  only  once. 

4.  Security/Safety 

Finally  there  is  the  security  advantage.  As  long  as  testing  is  being  conducted 
in  a  laboratory,  external  security  vulnerability  is  reduced.  Range  safety  costs 
and  concerns  that  occur  in  open-air  testing  are  also  eliminated. 
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To  accomplish  valid  and  accurate  laboratory  testing,  careful  attention  must  be 
paid  to  the  test  configuration.  For  HITL  tests  it  must  be  ensured  that  the  laboratory 
equipment  can  recreate  the  dynamics  and  conditions  that  will  be  experienced  in  the 
real  world.  Often  it  is  easy  to  switch  between  configurations  involving  hardware  and 
software.  This  requires  accurate  software  models  of  the  hardware  being  simulated. 
This  is  where  GLFF  networks  discussed  in  this  work  enters  the  application  front. 

To  gain  a  better  understanding  of  the  systems  the  GLFF  networks  must  model,  a 
more  detailed  breakdown  of  a  typical  missile  system  is  now  given.  The  functional 
components  of  the  missile  are  discussed  as  well  as  how  they  are  integrated  into  a 
HTTL  test  configuration.  The  GLFF  networks  can  then  be  used  to  approximate  some 
of  these  components. 

5.2.2  Functional  Breakdown  of  a  Missile  System 

This  section  provides  a  discussion  of  components  (subsystems)  most  often  found 
in  missile  systems.  Figure  5-7  below  illustrates  a  high-level  functional  break  down  of 
a  typical  missile  system  [Gar80]. 
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Figure  5-7:  Missile  functional  block  diagram. 


93 

Seeker/Guidance  Modules.  It  can  be  said  that  the  seeker  is  the  eye  and  the 
guidance  module  is  the  brain  of  a  weapon  system.  It  is  here  the  missile  trajectory  is 
commanded  and  controlled.  For  missiles  that  are  considered  lockon  after  launch,  the 
guidance  module  may  begin  in  a  target  search  mode.  In  this  case  it  drives  the  seeker 
to  search  for  a  target  in  some  pre-programmed  pattern.  During  this  phase  commands 
may  be  sent  to  the  autopilot  to  have  it  execute  a  midcourse  maneuver  such  as  pitching 
the  weapon  up  in  order  to  gain  altitude  for  target  search  and  future  end  game 
engagement.  If  target  energy  is  found  the  guidance  module  may  transition  into  an 
acquisition  mode  to  narrow  the  search  field  of  the  seeker  and  attempt  to  select  a  single 
target.  At  this  point  it  begins  target  track  and  will  try  to  keep  the  seeker  pointed  in  the 
direction  of  the  target.  During  track  mode,  guidance  (acceleration)  commands  are 
sent  to  the  autopilot  in  accordance  with  the  guidance  law  being  executed  (proportional 
or  pursuit  navigation).  All  the  while  this  module  performs  seeker  head  stabilization 
and  directional  pointing  to  maintain  a  continuous  track  of  the  target  until  impact. 

Autopilot.  The  job  of  the  autopilot  is  to  maintain  a  stable  flight  and  deliver  on  the 
commands  received  from  the  guidance  module.  To  perform  its  task,  the  autopilot 
accepts  input  from  various  sources  and  sensors.  The  three  most  common  are: 
guidance  module,  accelerometers,  gyros.  The  guidance  module  tells  the  autopilot  that 
a  certain  amount  of  acceleration  is  required  to  maintain  the  desired  trajectory.  Upon 
receiving  a  command  for  a  specific  acceleration,  the  autopilot  issues  commands  to  the 
actuator  to  drive  the  control  surfaces  (fins  or  canards)  in  the  proper  direction. 
Measurements  are  taken  from  accelerometers  to  determine  whether  the  weapon  has 
achieved  the  desired  acceleration.    Gyro  rotational  rate  and  angle  measurements  are 
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taken  to  ensure  proper  damping  and  to  keep  the  motion  of  the  missile  in  a  stable 
condition. 

Actuator  and  Control  Surfaces.  The  purpose  of  the  actuator  is  to  accept 
commands  from  the  autopilot  and  in  turn  drive  the  control  surfaces  to  achieve  the 
dynamics  needed.  As  the  fins  are  controlled  into  the  wind  the  missile  body  will  rotate 
and  accelerate  hopefully  in  a  stable  manner  sufficient  to  maintain  a  target  track. 
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Figure  5-8:  Closed-loop  missile  simulation. 


5.2.3  FUTL  Closed-Loop  Missile  Simulations 

In  this  section  FUTL  closed-loop  missile  simulations  will  be  discussed.  Figure  5-8 
above  illustrates  a  typical  test  configuration.  Usually  the  seeker,  guidance  module, 
and  autopilot  are  hardware  components  and  the  functionality  of  these  systems  is  the 
subject  of  the  test.  Special  laboratory  equipment  is  used  to  generate  the  target  scene. 
This  is  dependent  upon  the  kind  of  sensor  in  the  seeker  head.  The  remaining  blocks 
are  simulated  in  software.    For  a  FOTL  closed-loop  simulation,  the  seeker,  guidance 
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module,  and  autopilot  are  mounted  on  a  flight  motion  simulator  (FMS).  This  FMS  is 
commanded  to  recreate  the  body  dynamics  that  the  weapon  would  usually  experience 
in  a  real  flight.  Clearly  the  ability  of  the  FMS  to  create  accurate  body  motion  is 
critical  to  a  valid  FflTL  simulation.  To  understand  the  performance  capability  of  the 
FMS  an  accurate  software  model  of  this  system  is  desired.  This  model  can  be  used  to 
study  the  limitations  of  the  FMS  without  subjecting  the  hardware  to  harsh 
experimentation.  Also,  the  all-digital  simulation  can  be  modified  to  provide  for  body 
rate  updates  either  directly  computed  or  filtered  by  the  FMS.  In  this  way,  effects  of 
the  FMS  on  the  accuracy  of  the  weapon  simulation  can  be  determined.  A  discussion 
is  now  given  of  the  FMS  system  that  will  be  modeled.  Its  functional  composition  is 
given  along  with  its  expected  performance  specification. 

5.3  Modeling  the  Step  Response  of  a  Carco  5-Axis  FMS 

The  flight  motion  simulator  (FMS)  under  study  is  a  CARCO  model  S-458R-5E 
five-axis  system.  This  simulator  is  a  precision,  electro-hydraulic  positioner  consisting 
of  three  basic  elements:  a  hydraulically  driven  flight  motion  table,  a  control  console 
with  interface  electronics,  and  a  hydraulic  power  supply. 

The  FMS  structure  provides  five  axes  of  rotational  motion,  using  precision 
electro-hydraulic  actuators.  The  inner,  middle,  and  outer  gimbals  of  the  three-axis 
missile-motion  table  are  designed  to  roll,  pitch,  and  yaw  the  weapon.  The  two  large 
outer  gimbals  are  used  as  target  axes  and  are  designed  to  provide  azimuth  and 
elevation  maneuvers  for  the  target  system. 


96 

The  electronics  console  houses  the  electronic  servo  amplifiers,  power  supplies  and 
other  control  electronics  needed  to  regulate  the  table  motion.  It  provides  for  both 
analog  and  digital  closed-loop  control.  A  high-speed  parallel  interface  allows  remote 
control  and  monitor  of  the  FMS  from  a  digital  simulation  computer.  This  digital 
interface  allows  for  real-time  DMA  transfers  at  a  rate  of  1  frame  per  millisecond  (1 
KHz  digital  update). 

The  hydraulic  unit  provides  the  power  needed  for  the  high  dynamic  performance 
of  the  five-axis  table.  It  consists  of  a  motor,  pump  and  reservoir. 

As  delivered,  the  FMS  system  has  the  following  capability  using  a  150-lb  load  16" 
in  diameter  30"  long  mounted  on  the  roll  plate.  The  FMS  is  illustrated  in  Figure  5-9. 


Figure  5-9:  Carco  5-axis  FMS. 
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The  performance  specifications  for  each  gimbal  axis  are  listed  in  Table  5-4.   It  is 
this  closed-loop  system  that  will  be  modeled  with  several  GLFF  networks. 

Table  5-4:  Carco  FMS  performance  specifications 


Roll 
Axis 

Yaw 
Axis 

Pitch 
Axis 

Azimuth 
Axis 

Elevation 
Axis 

Max  position  (deg) 

Continuous 

±45 

±60 

±45 

±45 

Max  velocity  (deg/sec) 

1,000 

300 

300 

100 

10 

Max  acceleration  (deg/sec  ) 

20,000 

15,000 

15,000 

12,000 

12,000 

A  typical  closed  loop  system  (plant  plus  feedback  controller)  is  shown  in  Figure 
5-10.  When  identification  of  only  the  plant  is  of  interest,  model  design  could  be 
carried  out  with  an  open-loop  experiment.  If  the  controller  design  is  known,  it  is 
possible  to  model  the  plant  in  a  closed-loop  experiment  [Wan94a].  This  is  often  more 
desirable  than  open-loop  system  identification  because  it  causes  less  disruption  to  the 
operation  of  the  system  since  disabling  the  feedback  is  not  necessary.  If  the  open- 
loop  system  is  unstable  or  poorly  damped  then  open-loop  experiments  can  not  be 
performed  without  concern  for  safety  and  impending  system  damage. 
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Figure  5-10:  Typical  closed-loop  system. 

In  this  thesis,  it  is  a  model  of  the  entire  closed-loop  system  that  is  of  interest.  This 
is  often  needed  when  the  system  is  to  be  treated  as  a  sub-component  of  a  larger 
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configuration.  In  our  case,  the  system  is  a  hardware  device  that  is  to  be  driven  by  a 
digital  computer.  Since  the  system  is  expensive,  prolonged  and  abusive  use  is  highly 
undesirable.  If  it  is  possible  to  model  the  system  in  software,  then  preliminary 
development  and  debugging  of  the  computer  algorithms  and  hardware  interface  can 
be  performed  using  a  digital  system  model.  Once  the  software  has  been  tested,  the 
hardware  system  can  then  replace  the  digital  model.  In  this  way,  hardware  utilization 
is  kept  to  a  minimum.  Hence,  an  accurate  model  of  the  closed-loop  system  is  desired. 

5.3.1  Discrete-Time  GLFF  Models  of  the  FMS 

Figure  5-11  illustrates  the  general  form  of  a  discrete-time  GLFF  network 
realization.  Three  transfer  function  kernels  (Gamma,  Laguerre,  Legendre)  will  be 
considered  for  this  application. 


Figure  5-11:  Discrete-time  GLFF  network  block  diagram. 


The  GLFF  formulation  of  the  discrete-time  Gamma,  Laguerre,  and  Legendre 
kernels  were  given  in  Table  3-1.  These  kernels  result  in  the  network  realizations 
illustrated  in  Figures  5-12,  5-13,  and  5-14.  Note  that  in  all  cases  //01(z,/l)  =  1. 
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Figure  5-12:  Gamma  GLFF  Form. 
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Figure  5-13:  Laguerre  GLFF  Form. 
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Figure  5-14:  Legendre  GLFF  Form. 
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5.3.2  Step  Response  Modeling  Procedure 

As  stated  above,  the  inner  three  FMS  axes  carry  a  weapon  under  test  and  induce 
body  pitch,  yaw,  and  roll  motion.  A  target  may  be  presented  on  the  outer  azimuth  and 
elevation  axes.  Each  axis  is  independently  controlled  to  produce  the  proper  dynamic 
behavior  that  the  weapon  would  sense  if  it  were  flying  in  the  real  world.  Obviously 
there  is  much  work  required  in  developing  a  valid  HITL  weapon  simulation.  During 
the  HITL  simulation  development  cycle  a  software  model  of  the  FMS  can  be  used  in 
the  place  of  the  real  system.  This  will  eliminate  unnecessary  wear-and-tear  on  the 
hardware.  Hence,  a  good  digital  model  of  the  FMS  is  the  goal.  A  procedure  is  now 
described  that  will  employ  GLFF  networks  to  accomplish  this  task.  The  model 
parameters  are  optimized  using  measured  step  response  data.  This  data  was  produced 
by  subjecting  the  FMS  to  an  analog  step  input.  The  input  signal  and  FMS  response 
were  digitized  and  collected  at  a  1  KHz  sample  rate.  The  following  analysis  will 
concentrate  on  the  pitch  axis  only.  The  same  procedure  was  applied  to  the  other  axes. 

In  addition  to  the  step  response  data,  FMS  pitch,  yaw,  and  roll  gymbal  commands 
and  responses  were  collected  during  a  typical  missile  flyout  engagement.  The  gymbal 
commands  were  also  fed  into  several  GLFF  models.  The  model  responses  were 
compared  to  the  FMS  signals  to  characterize  the  accuracy  with  which  the  GLFF 
networks  can  approximate  the  FMS  under  normal  operations.  This  engagement  test 
was  conducted  using  an  AGM-65G  Maverick  missile.  The  AGM-65G  is  a  short 
standoff-range  launch-and-leave  air-to-ground  missile.  It  is  a  roll-stabilized, 
cruciform-wing,  tail-controlled  system  with  an  imaging  infrared  sensor.  The  scenario 
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chosen  launches  the  missile  at  approximately  200  meters  altitude  and  2000  meters 
down  range  from  a  stationary  ground  target.  The  results  of  this  test  are  presented  with 
the  step  response  data  in  section  5.3.3. 

The  approach  taken  to  optimize  a  discrete-time  GLFF  network  parameter  set  using 
step  response  data  of  the  FMS  system,  H,  is  illustrated  in  Figure  5-15.  Since  the  goal 
is  a  network  that  models  the  impulse  response  of  H  via  step  response  data  a  derivative 
is  necessary;  hence,  precautions  must  be  taken  to  ensure  the  step  response,  g(t),  has 
low  signal-to-noise  ratio.  Another  approach  is  to  transfer  the  derivative  from  the 
response  signal  to  the  network  kernels.  This  and  other  methods  were  discussed  in 
Chapter  3.5. 


es(t) 


Figure  5-15:  Step  response  modeling  structure. 

Because  of  the  derivative  of  the  step  response  data,  the  parameter  adaptation 
process  must  work  to  reduce  the  overall  noise  effects.  The  process  involves  selecting 
the  parameter  set  p  =  {N,w,X\  in  a  two-stage  fashion. 

Step  1:  Fix  the  order  N  and  determine  the  time  scale  parameter,  A,  and  weight  vector, 
vv,  that  minimizes  the  integrated  squared  error,  J,  of  the  output  of  the  system  and 
model  impulse  responses. 
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Step  2:  Let  A  adapt  to  minimize  the  squared  error  of  the  output  of  the  system  and 
model  step  responses.  In  this  way,  the  optimal  weight  vector  represents  the 
impulse  response  while  the  optimal  time  scale  allows  for  a  model  match  of  the 
step  response.  In  a  noise  free  environment,  A  would  in  fact  be  the  optimal  time 
scale  minimizing  the  impulse  response  squared  error. 


This  method  is  now  demonstrated.  The  test  case  illustrates  modeling  the  FMS 
pitch  axis.  The  set  point  (step  input)  is  1  volt  (6  Degrees).  The  step  input  and  pitch 
axis  step  response  are  shown  in  Figure  5-16. 
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Figure  5-16:  Pitch  axis  step  response. 
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5.3.3  Results  (FMS  Step  Response  Modeling) 

GLFF  model  parameters  were  optimized  using  the  procedure  described  above  for 
the  first  eight  orders  of  Gamma,  Laguerre,  and  Legendre  networks.  Tables  5-6,  5-7, 
and  5-8  summarize  the  numerical  results.    Here,  A,     is  the  optimal  time  scale  from 

step  1.  J Inun  is  the  minimum  squared  error  resulting  from  the  Wiener-Hopf  solution 
of  the  impulse  responses,  Jlmin  =  J(A,  ).  Likewise,  As  t  is  the  optimal  time  scale 
from  step  2.  JSnin  is  the  minimum  squared  error  resulting  from  the  adaptation  of  the 
step  responses,  J Smm  =  J(AS  ,).    As  can  be  seen  from  the  tables,  all  three  kernels 

yield  similar  capabilities.  Because  of  the  completeness  of  the  GLFF  structure, 
increasing  the  order  of  the  network  will  always  produce  either  the  same  or  smaller 

Imin  ' 

For  each  (Gamma,  Laguerre,  Legendre)  network  squared  error  performance 
functions  (Figures  5-21,  5-26,  and  5-31),  step  responses  (Figures  5-22,  5-27,  and  5- 
32),  and  impulse  responses  (Figures  5-23,  5-28,  and  5-33)  are  plotted  for  models  of 
order  1  to  8.  For  the  8th  order  models,  the  step  response  and  associated  error  (Figures 
5-24,  5-29,  and  5-34),  and  flyout  response  and  associated  error  (Figures  5-25,  5-30, 
and  5-35)  are  plotted.  Finally  the  8l  order  step  responses  are  scaled  to  degrees  and 
plotted  in  Figure  5-36  along  with  the  scaled  errors  in  Figure  5-37.  Notice  that  the 
maximum  error  is  no  more  than  0.15°  for  a  6°  step  (less  than  3  percent). 

Once  the  time  scale  was  found  that  produced  the  minimum  squared  error  for  the 
impulse  response,  //„„■„,  it  was  allowed  to  adapt  to  search  for  the  time  scale  yielding 
the  minimum  squared  error  for  the  step  response.   These  are  slightly  different  due  to 
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the  noise  in  the  impulse  response.   Choosing  JSmm  -  -30dB  as  the  error  specification 


i  tii 


,ld 


rd 


a  4     order  Gamma,  3     order  Laguerre,  and  3     order  Legendre  model  would  be 


chosen. 


Figures  5-38  and  5-39  plot  the  transfer  function  magnitude  and  phase  of  the 
system  and  models.  They  match  very  well  out  to  30  Hz.  The  pitch  axis  frequency 
response  under  the  specified  load  condition  is  20  Hz.  The  system  response  increase 
from  60  -  100  Hz  is  a  result  of  this  being  a  closed-loop  system. 

To  provide  a  comparison  of  the  GLFF  networks  with  the  conventional  tap-delay 
line,  Table  5-5  lists  the  results  of  several  FIR  model  performance  parameters.  Figures 
5-17  through  5-20  plot  the  step  response,  impulse  response,  step  response  error,  and 
flyout  response  using  a  201 -tap  delay  line.  A  large  order  FIR  model  is  required  to 
achieve  good  results.  This  is  a  disadvantage  due  to  increased  computational 
requirement  as  well  as  the  system  delay  if  the  model  had  to  operate  in  a  real-time 
environment. 


Table  5-5:  FIR  optimization  characteristics. 


Order 

Imin 

J, 

Snun 

51 

-3.2 

-10.4 

101 

-12.3 

-13.1 

151 

-14.4 

-19.4 

201 

-14.9 

-21.9 
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Figure  5-17:  FIR  network  impulse  response  (201  taps). 
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Figure  5-18:  FIR  network  step  response  (201  taps). 
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FIR  Model  Pitch  Axis  Step  Response  (order=201 ) 
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Figure  5-19:  Step  response  and  associated  error  for  201-tap  FIR  network. 
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Figure  5-20:  Flyout  response  and  associated  error  for  201-tap  FIR  network. 
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DT  Gamma  Model  Pitch  Axis  Performance  Functions  (orders=1  to  8) 
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Figure  5-21:  Gamma  network  performance  functions. 


Table  5-6:  Gamma  network  optimization  characteristics. 


Order 

K> 

Imin 

^Sopt 

J, 

Smin 

1 

0.0135 

-2.9 

0.0205 

-19.6 

2 

0.0380 

-5.9 

0.0530 

-22.5 

3 

0.0620 

-7.6 

0.0260 

-29.9 

4 

0.0415 

-9.4 

0.0420 

-33.0 

5 

0.0590 

-10.8 

0.0405 

-33.0 

6 

0.0715 

-11.1 

0.0295 

-36.0 

7 

0.0895 

-11.7 

0.0415 

-40.7 

8 

0.1050 

-11.9 

0.0500 

-41.7 

108 


1.4 
1.2 

DT  Gamma  Model  Pitch  Axis  Step  Response  (orders=1  to  8) 

i             i             i             i             i             i             i             i 

1 

■     y/^^ 

gGAM(t)  (volts) 
o          o 

CD               CO 

I               : 

S  0.4 

71 

0.2 

'I 

0 

i      i      i      i      i      i      i      i 

0.1  0.2         0.3         0.4         0.5         0.6  0.7         0.8         0.9  1 

t  (sec) 

Figure  5-22:  Gamma  network  step  responses. 
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Figure  5-23:  Gamma  network  impulse  responses. 
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DT  Gamma  Model  Pitch  Axis  Step  Response  (order=8) 
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Figure  5-24:  Gamma  order  8  step  response  and  associated  error. 
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Figure  5-25:  Gamma  order  8  flyout  response  and  associated  error. 
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DT  Laguerre  Model  Pitch  Axis  Performance  Functions  (orders=1  to  8) 
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Figure  5-26:  Laguerre  network  performance  functions. 


Table  5-7:  Laguerre  network  optimization  characteristics. 


Order 

^lopt 

Imin 

^■Sopl 

1 

0.9625 

-5.8 

0.9475 

-22.5 

2 

0.9390 

-7.5 

0.9745 

-29.9 

3 

0.9595 

-9.3 

0.9585 

-32.9 

4 

0.9415 

-10.8 

0.9605 

-33.0 

5 

0.9290 

-11.1 

0.9710 

-35.7 

6 

0.9120 

-11.6 

0.9585 

-40.3 

7 

0.8950 

-11.9 

0.9500 

-41.6 

8 

0.9375 

-12.1 

0.9505 

-41.7 
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Figure  5-27:  Laguerre  network  step  responses. 
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Figure  5-28:  Laguerre  network  impulse  responses. 
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Figure  5-29:  Laguerre  order  8  step  response  and  associated  error. 
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Figure  5-30:  Laguerre  order  8  flyout  response  and  associated  error. 
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DT  Legendre  Model  Pitch  Axis  Performance  Functions  (orders=1  to  8) 


Figure  5-31:  Legendre  network  performance  functions. 


Table  5-8:  Legendre  network  optimization  characteristics. 


Order 

^lopl 

Imin 

2 

"'Sopi 

Smin 

1 

0.9729 

-5.6 

0.9615 

-22.2 

2 

0.9642 

-6.9 

0.9863 

-29.2 

3 

0.9822 

-8.8 

0.9812 

-33.1 

4 

0.9775 

-10.3 

0.9770 

-34.3 

5 

0.9846 

-10.6 

0.9819 

-35.1 

6 

0.9891 

-11.0 

0.9883 

-38.8 

7 

0.9690 

-11.4 

0.9866 

-40.8 

8 

0.9854 

-11.8 

0.9848 

-42.2 
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Figure  5-32:  Legendre  network  step  responses. 
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Figure  5-33:  Legendre  network  impulse  responses. 
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DT  Legendre  Model  Pitch  Axis  Step  Response  (order=8) 
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Figure  5-34:  Legendre  order  8  step  response  and  associated  error. 
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Figure  5-35:  Legendre  order  8  flyout  response  and  associated  error. 
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Figure  5-36:  FMS  and  8th  order  Gamma,  Laguerre,  Legendre  network  step  responses. 
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Figure  5-37:  8    order  Gamma,  Laguerre,  Legendre  networks  step  response  errors 
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Laguerre(8)  /  Gamma(8)  /  Legendre(8)  Transfer  Function  Estimates 
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Figure  5-38:  Transfer  function  magnitude  estimates  for  FMS  and  8l  order  networks. 
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Figure  5-39:  Transfer  function  phase  estimates  for  FMS  and  8l  order  networks. 
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5.3.4  Continuous-Time  GLFF  Models  of  the  FMS 


In  this  section  the  GLFF  modeling  capability  will  be  further  illustrated  by 
constructing  a  continuous-time  Laguerre  model  of  the  FMS  using  the  sampled  step 
response  data.  The  general  form  of  a  continuous-time  Laguerre  network  using  the 
Laguerre  kernel  is  given  in  Figure  5-40. 
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Figure  5-40:  Continuous-time  GLFF  network  block  diagram. 

The  optimization  procedure  employed  is  that  which  was  derived  in  Chapter  4. 
Namely,  the  analytic  performance  function,  Ja,  will  be  computed  and  the  roots,  Ar,  of 
2Im[7a]  located.  A  zero  crossing  detection  algorithm  was  developed  to  numerically 
evaluate  2Im[7fl]  for  these  roots.  The  squared  error  is  calculated  for  each  root  and  the 
time  scale  producing  the  minimum  is  selected  as  the  nearest  ZCa. 

A^=argfmin{2Re[7fl(2r)]}j  (5.5) 

Figure  5-41  depicts  the  squared  error  performance  functions  for  Laguerre  networks 
from  order  0  to  8.    For  each  order  network,  Table  5-9  provides  a  comparison  of  the 


optimal  time  scale,  Aopl,  the  nearest  ZCa  time  scale,  /lnZca,  and  their  respective  squared 


119 

error  values.  Figures  5-42  and  5-43  show  the  squared  error  and  imaginary  component 
of  the  analytic  performance  function.  Recall  J  is  embedded  within  Ja  and  is 
calculated  by 

J  =  (daeR  ~  2wlRWaR)ldaeR  (5-6) 

Finally,  Figure  5-44  graphs  the  step  responses  and  step  response  errors  for  8th 
order  models  using  Aop,  (light  curve)  and  A„zca  (dark  curve)  as  the  time  scale.  Figure 
5-45  graphs  the  responses  and  associated  errors  of  these  models  when  driven  by  the 
AGM-65G  flyout  commands. 
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Figure  5-41:  Laguerre  network  performance  functions. 
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Table  5-9:  Laguerre  network  optimization  characteristics. 


Order 

A  apt 

/lnZCa 

JUopt) 

JttnZCa) 

0 

13.5 

11.0 

-2.6 

-2.6 

1 

38.5 

41.0 

-5.3 

-5.2 

2 

63.0 

69.0 

-6.7 

-6.5 

3 

42.0 

39.0 

-8.0 

-7.9 

4 

60.5 

60.0 

-9.0 

-9.0 

5 

74.0 

55.0 

-9.2 

-9.1 

6 

92.5 

39.5 

-9.5 

-9.3 

7 

111.0 

49.0 

-9.6 

-9.5 

8 

65.0 

63.0 

-9.8 

-9.7 
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Figure  5-42:  J  and  2Im[7a]  for  FMS  pitch  axis  (orders  =  0  to  3). 
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Figure  5-43:  J  and  2Im[/a]  for  FMS  pitch  axis  (orders  =  4  to  7). 
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5.3.5  Conclusion  (Step  Reponse  Modeling) 

The  purpose  of  this  chapter  has  been  to  review  and  illustrate  the  applicability  of 
GLFF  networks  in  the  area  of  system  modeling.  Focus  was  on  the  following 
concepts: 

1.  A  methodology  for  determining  GLFF  network  realizations  of  the  impulse 
response  of  a  closed-loop  system  using  measured  step  response  data  is 
possible.  The  procedure  is  simple  and  can  be  practically  implemented  in  a 
laboratory  setting. 

2.  A  Carco  five-axis  flight  motion  simulator,  used  in  weapon  system  hardware  - 
in-the-loop  test  and  evaluation,  can  be  modeled  very  accurately  using  low 
order  GLFF  networks.  This  can  be  quite  useful  since  a  digital  model  of  the 
FMS  could  be  substituted  for  the  hardware  when  missile  flyout  models  are 
under  development,  or  when  FMS  response  issues  are  raised  and  the  hardware 
is  unavailable. 

5.4  Other  Applications 

When  conducting  laboratory  HITL  closed-loop  simulations,  some  of  the  weapon's 
subsystems  such  as  the  actuators,  gyros,  and  accelerometers  must  be  modeled.  If 
hardware  actuators  are  available,  they  must  be  pressurized  and  loaded  to  provide 
accurate  fin  position.  This  becomes  impractical  to  do  on  a  simulation-by-simulation 
basis.  If  a  single  step  response  test  was  performed  and  the  tools  described  in  this 
research  employed,  then  a  GLFF  network  could  be  constructed  to  accurately  represent 
and  replace  the  actuators.  Gyros  are  also  difficult  to  work  with  in  HITL  testing  and 
they  are  usually  modeled.  Accelerometers  can  not  be  exercised  because  no  body 
translational  motion  is  created  in  the  laboratory.  Again,  if  hardware  is  available  and 
can  be  statically  evaluated  to  provide  step  response  data,  GLFF  networks  could  be 
constructed  to  model  these  subsystems.  This  offers  an  advantage  over  classical 
models  because  the  GLFF  networks  provide  a  representation  of  the  real  systems  as 


124 

compared  to  the  theoretical  approach  where  the  internal  parameters  (gains,  time 
constants,  damping  coefficients,  etc.)  are  considered  to  be  known  and  assumed  to  be 
operating  with  exact  precision. 

5.5  Conclusion  (Applications) 

This  chapter  has  discussed  various  applications  of  the  GLFF  model  class.  A 
historical  overview  of  the  many  kinds  of  problems  that  have  been  solved  was  first 
presented.  Next,  several  examples  were  given  to  illustrate  how  these  models  can  be 
used  to  create  networks  that  approximate  systems  when  only  step  response  data  is 
available.  Having  the  model  parameters  the  impulse  and  step  responses  were 
recreated  and  compared  to  the  true  systems. 

Since  the  "real  world"  applications  of  interest  are  systems  and  subsystems 
involved  in  test  and  evaluation  of  guided  weapons  the  T&E  process  was  described. 
This  included  the  importance  of  laboratory  tests  such  as  all-digital  and  HITL 
simulations.  A  HITL  test  configuration  was  given  and  was  used  to  identify  various 
systems  that  are  good  candidates  for  GLFF  networks.  Primarily,  the  FMS  system  was 
studied  and  several  GLFF  networks  were  created  using  sampled  step  response  data. 


APPENDIX  A 
ORTHOGONAL  SIGNAL  SETS 


A  set  of  complex-valued  functions  of  the  real  variable  t  (A  is  considered  fixed 
here),  {0n(t,A)}~   ,  «  =  0,1,2,  ...  is  called  orthogonal  [Su71]  over  the  range  (a,b)  if 

\</>m(t,A)fn(t,A)dt  =  \  (A.l) 

J  nonzero     m  =  n 

a  L 

The  set  is  called  orthonormal  over  the  range  (a,b)  if 

U  (U^ca^H,  (A.2) 

J  1     m  =  « 
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APPENDIX  B 
HILBERT  SPACE  CONCEPTS 


Definition:  (Inner  Product) 

Let  E  be  a  complex  vector  space.  A  mapping 

(;»):ExE->€  (B.I) 

is  called  an  inner  product  in  E  if  for  any  x,y,z  £  E  and  any  a,(3  6  C  the  following 
conditions  are  satisfied: 

A.  (x,y)  =  (y,x}* 

B.  (ax  +  0y,z)  =  a{x,z)  +  0(y,z)  (B.2) 

C.  (x,  x)  >  0  and  (jc,  x)  =  0  implies  x  -  0 

Definition:  (Inner  Product  Space) 

A  vector  space  with  an  inner  product  is  called  an  inner  product  space. 

Definition:  (Norm  in  an  Inner  Product  Space) 

The  norm  in  an  inner  product  space  E  is  defined  as  the  functional 
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Remark: 

The  space  Lr([a,b])  of  all  Lebesgue  square  integrable  functions  on  the  interval  [a,b] 

with  the  inner  product 

b 

(f,g)  =  \f(t)g(t)dt  (B.4) 

a 

is  an  inner  product  space. 

Definition:  (Cauchy  Sequence) 

A  sequence  [xn},  n=l,2,...,  is  a  Cauchy  sequence  if  for  every  £ >  0  there  exists  a 

number  M  such  that  be  -  x  II  <  £  for  all  n,  m  >  M. 


m  n\ 


Definition:  (Complete  Inner  Product  Space) 

An  inner  product  space  is  complete  if  every  Cauchy  sequence  from  the  space 

converges  to  an  element  in  the  space. 

Definition:  (Hilbert  Space) 

A  Hilbert  Space  is  a  complete  inner  product  space. 


APPENDIX  C 
S-DOMAIN  DERIVATION  OF  MINIMUM  J 


Using   Parseval's   theorem   we   can   represent   the   squared   error   performance 
function  in  the  s-domain  as 

j=  je\t)dt  =  ^§E(s)E(-s)ds  (C.l) 

—  c 

where  E(s)  =  D(s)-  Y(s)  and  C  is  a  contour  taken  in  a  CCW  direction  enclosing  the 
entire  left  hand  plane  (LHP)  in  s.  If  Y(s)  is  the  OGLFF  model  of  the  impulse  response 
then 

E(s)  =  D(s)-fjwn<&n(s,A)  =  D(s)-wTd>(s,A) 

n=0 

where  (C.2) 

wT  =  (w0,wl,...,wN)<md®(s,A)  =  (<S>0(s,A),<t>l(s,A),...,®N(s,A)jr 

So  the  performance  function  becomes 

J=^j[D(s)-wT^(s,A)][D(-s)-wT^(-s,A)}is 
c 

=  ^JD(s)D(-s)ds-wT^-JD(sm-s,A)ds  (C.3) 

c  c 

-wT -±-i D(-s)$>(s,A)ds  +  wT ^-§<&(s,A)®T (s,A)ds  w 
c  c 

Since  {0„(.s,A)}     is  an  orthonormal  sequence 

if$fl(a)^(-a)^=^  (c.4) 
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or 


^§<S>(s,A)®T(-s,A)ds=I  (C.5) 

c 

Also,  for  the  contour  C  we  have 

^§D(sm-s,A)ds  =  ^§D(-s)®(s,A.)ds  (C.6) 

c  c 

With  the  definition   de  =  \d'{t)dt  =  j^§D(s)D(-s)ds  which  is  the  energy  of  the 

0  c 

signal  d(t)  then 

J  =  de-wT^§D(s)<S>(-s,A)ds  +  wTw  (C.7) 

c 

For  a  fixed  A,  to  minimize  7  we  must  satisfy  Vw7  =  0  .  Hence 

c 
which  yields 

w  =  ^lD(s)®(-s,A)ds  (C.9) 

c 
as  the  optimal  weight  vector  that  minimizes  the  squared  error  performance  function. 
Substituting  w  into  the  equation  for  J  gives 

Jop,=de-wTw  (CIO) 


APPENDIX  D 
SQUARED  ERROR  STATIONARITY  USING  THE  LAGUERRE  MODEL 


The  necessary  conditions  for  stationarity  of  the  squared  error  performance 
function  when  using  a  Laguerre  model  are  now  derived  [Wan94b].  The  derivation 
makes  use  of  the  following  lemma,  easily  proven  using  induction. 


Lemma:  For  the  Laguerre  kernels  Ln(s,A)  =  4lA -,     A>0,n  =  0,1,... 

(s  +  A)"+ 

the  following  relation  is  true 
"  Ln{s,A)={-^^L^(s,A)-^-Ln_x{s,A\     L_,(s,A)  =  0  (D.l) 


dA    n  2k       "+1V  '    '     2/1 

In  Chapter  4.3.1  the  minimum  squared  error  performance  function  for  a  GLFF 
model  (using  the  Laguerre  kernel)  was  derived 

N 

J  =  de-wTw  =  de-^w;  (D.2) 


«=0 


where 


de=^JD(s)D(-s)ds 
c 

w  =  ^$D(s)L(-s,A)dS 
c 

w  =  (w0,wl,...,wN)T 

L(s,A)  =  (LX)(s,A),Ll(s,A),...,LN(s,A))T 


(D.3) 
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J  achieves  a  minimum  when  2^w^is  a  maximum.   So  the  stationary  points  are  found 


n=0 


by 


Af^(i)=2fw>a)^w=0 


^  „=0 


«=0 


(D.4) 


■B 


Now 


^-^^A* 


(D.5) 


Using  the  above  Lemma  and  simplifying  gives 
dwn{A)      )i  +  [ 


w„+1(A) w„.(A),     w,(A)  =  0 

<?/l  2A     "+1  2X    ""'  "'      J 


(D.6) 


So  V  w^  is  a  maximum  when 


n=0 


&  A-  n=0  n=0 


H  +  l  ,  3  v  "  f  3  v 

2A  2A    ""' 


=  0 


(D.7) 


or 


(^  +  l)w„U)WA,+ia)  =  0 


(D.8) 


Hence  the  optimal  value  of  X  is  either  a  root  of  wN (A )  =  0  or  a  root  of  wN+1  (A)  =  0. 
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